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Project Summary 


Over the past four years the Computational Mechanics Company, Inc. has been the principal 
investigator into the development of a new class of computational methods for modeling 
hypersonic viscous flows in two- and three-dimensional domains. The ultimate goal of this 
project is to provide NASA with a highly accurate computational tool for resolving fine 
scale flow features typical of shock wave interactions and strong viscous interaction regions 
in hypersonic flows. Toward this end, a research and development effort was put forth in 
the area of adaptive computational finite element methods for high speed flows based on 
unstructured mesh concepts and employing local estimates of the solution error to optimally 
change the computational grid to minimize the numerical error. The approach developed 
here, which combines /i-adaptive and anisotropic p-enrichment mesh modification procedures, 
has proven for certain classes of problems to provide exponential rates of convergence of the 
solution using possibly an order of magnitude fewer degrees of freedom than conventional 
methods while obtaining the same level of computational accuracy. 

During the first phase of this project (years 1-3), work was focused on a number of 
research topics which were crucial to the success of h-p adaptive finite element methods for 
hypersonic flow's. Many of these issues were resolved and the results are presented again in 
detail in the body of this report for completeness. Summarizing some of the more significant 
developments of this first phase of the effort: 

• The development of the first h-p finite element data structure for quadrilateral and 
hexahedral elements. 

• The development of an h-p adaptive package that includes 

— A local error estimation capability for driving the adaptive strategy. 

- A spectral enrichment and /i-refinement methodology to change the structure of 
the computational grid. 

• The formulation of a generalized methodology for handling nodal point constraints 
for higher order polynomials. (Note that this difficulty does not occur with either 
/i-refinement or p-enrichment individually.) 

• The development of an algorithm for manual anisotropic enrichment of both two- 
dimensional and three-dimensional elements. 

• The formulation and implementation of a version of a preconditioned block Jacobi- 
GMRES method for higher order spectral elements. 
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• The formulation and implementation of a one-step Taylor-Galerkin solution algorithm 
for h-p methods. 

• The formulation and implementation of an implicit two-step Taylor-Galerkin solution 
algorithm which solves first an Euler step followed by a viscous step. 

• An investigation of artificial dissipation mechanisms appropriate for h-p adaptive com- 
putational methods with possibly highly distorted elements or elements with high as- 
pect ratios. 

In addition to these efforts on research-oriented topics, a user friendly, graphics oriented, 
interface for the two- and three-dimensional codes was developed. This interface includes 
both a batch and interactive option and full graphics capabilities for displaying the solu- 
tions, extracted quantities, and plots of the local estimates of the computational error. A 
simple grid file interface was also developed to read in neutral grid files generated either 
by the GAMMA2D or GAMMA3D codes, developed in-house at COMCO, or by PATRAN. 
(The formats necessary for PATRAN and the neutral file interface are discussed in the grid 
generation section of the user manual.) 

The second phase of this effort, conducted over the last 12 months, has focused on two 
special research and development topics which are in general related to the performance of 
the flow solver. These topics include: 

• The development of implicit/explicit computational methods (in two dimensions) for 
integrating the Navier Stokes equations forward in time. 

• The development of computational methodologies and algorithms which provide for 
automated directional p-enrichment of the computational mesh. 

A third topic on which we dedicated considerable computer resources was a continued 
investigation of artificial dissipation mechanism for h-p adaptive computational methods. 
In particular, numerous test cases were run for the Mach 14 Holden problem using various 
artificial mechanism to determine an appropriate model for capturing the recirculation bubble 
and other fine features of the solution. 

All of the capabilities and options developed during the first and second phases of the 
effort have been implemented and/or tested in a two-dimensional and/or a three-dimensional 
finite element code. Using these codes a number of test cases have been solved to verify the 
functionality of the software. These benchmark cases include flow past a blunt body with 
an incident shock to produce a Type IV (Edney) interaction, a Holden problem with inflow 
Mach number 14, and a rearward facing step with a strong expansion. Experimental data 
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on each of these benchmark problems is available and comparisons with this data are made 
throughout the results section of this report. 
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1 Introduction 


The commitment to develop the National Aero-Space Plane and Maneuvering Reentry Ve- 
hicles has generated resurgent interest in the technology required to design structures for 
hypersonic flight. As these vehicles cruise, accelerate and/or decelerate in the atmosphere, 
highly complex patterns of shock wave interactions and shock wave boundary layer inter- 
actions develop which produce severe local pressures and extreme local heating rates. To 
provide adequate safety factors in the thermal-structural designs an accurate determination 
of the aerothermal loads, especially in the local areas of strong shock interactions and strong 
viscous interactions such as interacting shock waves on a leading edge, shock/boundary layer 
interactions, corner flows, etc., is required. In general, ground based test facilities can pro- 
vide only limited data for the expected flight conditions, at a considerable cost, and thus 
designers must depend heavily on analytically predicted aerothermal loads. The analyti- 
cal methodologies for these predictions must be accurate (to within some predetermined 
measure), robust, computationally economical, and geometry independent. 

For a designer, the prediction of aerothermal loads of hypersonic structures is one of the 
most challenging problems in computational fluid dynamics. It requires flow analysis and 
shock prediction at very high Mach numbers, as well as a realistic calculation of aerother- 
mal loads. These tasks required, in general, advanced computational strategies, extremely 
accurate discretization techniques, and powerful postprocessing capabilities. 

During the first two phases of this project we have developed, at the Computational 
Mechanics Company, an implicit/explicit, anisotropic h-p finite element methodology for 
the analysis of high speed flows and the prediction of aerothermal loads. The basic idea 
of the h-p version of the finite element method is to combine local mesh refinement (an 
A-method) with anisotropic polynomial enrichment (a p-method) in order to achieve con- 
vergence rates not attainable with fixed mesh methods or with any of the above methods 
applied separately. In practical terms, the h-p method means maximum numerical accuracy 
at a minimal computational cost. 

The remainder of this report summarizes the results of the first and second phases of a 
four year research and development effort oriented toward the resolution of several theoretical 
and computational issues related to the above problems. In Section 2 the basic formulation 
of the Navier-Stokes equations which governs the compressible viscous flow is presented. 
This is followed by a detailed discussion of a general family of implicit Taylor Galerkin algo- 
rithms for solving the compressible flow equations, and a review of the numerical boundary 
conditions associated with the algorithms. The final part of this section provides a brief 
discussion about artificial dissipation mechanisms which may be more appropriate for highly 
distorted elements or elements with high aspect ratios. The Taylor- Galerkin formulation 
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and the general boundary condition treatment presented in Section 2 is excerpted from our 
published work on these results (e.g. [8, 10, 31]). The next section, Section 3, provides a 
detailed review of the h-p finite element methodology. Included in this review is a discussion 
of the data structure, hierarchical shape functions, constrained approximation, and assump- 
tions and restrictions which have been incorporated into the h-p formulation. The work 
reported in Sections 3.3-3.11 is based on extensions of the general h-p adaptivity methods of 
Oden, Demkowicz, and Rachowicz, originally published in 1989 (see [9, 25]) and expanded 
and generalized to apply to thermomechanical load predictions in this project. Section 4 
follows with a discussion of the adaptive strategy which is used in conjunction with the 
mesh modification algorithms. Here the error estimates used to drive the adaptive package 
are outlined and the procedure for adapting the computational mesh is presented. Section 
5 follows with a detailed discussion of the implicit/explicit methodologies that have been 
implemented within the context of the two-dimensional code. In particular, in this section 
we review the basic concepts associated with using implicit/explicit methods and how one 
would select implicit end explicits zones within the computational domain. This section con- 
cludes with a brief overview of other computational procedures also required to efficiently 
implement an implicit/explicit methodology. 

The contents of this report includes a summary of the efforts completed during both the 
first and second phases of the development effort. During the second phase of the project, 
our efforts have focused on implicit/explicit methods, directional enrichment techniques and 
algorithms, and optimal artificial dissipation mechanisms. The details of this effort are 
provided in sections 2.2, 4.2, 5, and 7. The other basic sections are essentially the same as 
in previous reports except for various corrections. 

The next section, Section 6, presents some non-standard algorithms used in the element 
calculations and the postprocessing module. This is followed by Section 7, which presents 
the results of the numerical test cases that have been run over the course of this project. This 
includes both simple test problems used for verification purposes as well as the benchmark 
problems supplied by NASA-Langley. The final section provides a brief discussion of some 
possible directions for future research and development in the area of h-p adaptive methods 
for hypersonic flows. 
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2 Formulation of the Governing Equations and Nu- 
merical Algorithms 


This section provides a summary of the basic field equations and numerical algorithms that 
have been used to model the complex flow phenomena that is encountered in hypersonic 
flight conditions. We begin this section with a standard review of the basic field equations 
of compressible gas dynamics which are expressed in conservation form. These equations 
are then nondimensionalized in the usual fashion (see reference [1]) using the free stream 
density, velocity, Mach number, Reynolds number, etc. This procedure results in a system of 
evolution equations, for a set of nondimensionalized conservation variables, which assume es- 
sentially the same form as the original governing equations but with new material constants. 
The next two sections present the details of two algorithms for integrating the Navier-Stokes 
equations forward in time. The first algorithm is a one-step implicit/explicit method which 
is based on the general family of Taylor- Galerkin method with several parameters controlling 
the actual implicitness of the scheme. The second method is a two-step method which sepa- 
rates the Euler fluxes from the viscous contributions, resulting in an inviscid convection step 
and a viscous diffusion step. The next subsection summarizes the various classes of bound- 
ary conditions that arise from these two approaches and the final section discusses artificial 
dissipation mechanisms appropriate for h-p adaptive methods. (As a matter of notational 
simplicity all of the formulations and algorithms presented in this section are limited to the 
two-dimensional case.) 

2.1 Notation and Formulation of the Equations 

The compressible viscous flow of a calorically perfect gas is governed by the Navier-Stokes 
equations which may be conveniently written as 

U, t + Ff(U), i =F?(U), i (2.1) 

where U is the vector of conservation variables defined by 

U = {p,mi,m 2 ,e} T (2.2) 

and Ff{U) and F)' (U) are the inviscid and viscous fluxes, respectively. The indices t in (2.1) 
refer to the axis of a Cartesian coordinate system, a comma denotes partial differentiation, 
and the summation convention is applied. The Eulerian fluxes Ff (U) are defined as 
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Ff(U) 

flm 




mi 


+ P' 


raimj 


, (e + p) 


t)’ 


miirii m? , . m 2 V 

“T + p, (« + P -r) 
p p p / 


(2.3) 


Here p is the mass density, m, = pu, are the momentum components (with Ui the velocity 
vector components), e is the total energy per unit volume, and p is the thermodynamic 
pressure. Under the assumption that the fluid behaves as a perfect gas, the constitutive 
equation relating the pressure to the internal energy is given by 


p = (7 - l)i (2.4) 

where t is the internal energy per unit volume 

1 = e - \p( u \ + u l) ( 2 - 5 ) 

and 7 is the ratio of specific heats 


(2.6) 


c v 

In this expression C p is the specific heat at constant pressure and C v is the specific heat at 
constant volume. 

The viscous fluxes F ) ( U ) are defined by 


Fi (U) — (0, r iii r i2 1 ' r ij u j + 9i) ^ 7) 

F 2 [U) = (0, r 21 , r 22 , t 2j Uj + q 2 ) 

with the viscous stresses 77, related to the velocity gradient through the usual constitutive 
relation for a Newtonian fluid 


T tJ = p[u it j + «>,,•) 4- A Uk,kf>ij (2.8) 

and the heat flux is defined by Fourier’s law 

9, = xT'i (2.9) 

Here p and A are the dynamic viscosity and second viscosity coefficient, respectively, 6,7 is 
the Kronecker delta, T is the temperature, and « is the heat conduction coefficient. 
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In addition to the quantities defined above, we also introduce the Reynolds number Rei 
and the Prandtl number Pr 


Rti = 


pVL 


( 2 . 10 ) 


Pr = 


Cp/i 


( 2 . 11 ) 


K 

where L is the referential or characteristic length. These quantities will be used in the next 
section where a nondimensionalized form of the governing equations is discussed. 

As a final note, we will assume throughout this report that the Stokes relation is in effect 


3A -|- — 0 

and Southerland’s law relating the temperature to the dynamic viscosity holds 


j-3/2 

* = Cl T + C 2 

where Ci and C 2 are constants for a given gas. 


( 2 . 12 ) 


(2.13) 


2.1.1 Nondimensional Form of the Navier-Stokes Equations 

Following the procedure outlined by Anderson [1] we introduce the following scaling param- 
eters 

p ^ - free stream density 

- free stream velocity 
Too - free stream temperature 
Pea - free stream dynamic viscosity 
L - referential or characteristic length 

Using these scaling parameters, a number of nondimensional quantities may be defined 
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Voot 
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UL 

Ko 


P‘ 


P 

Hoc 


p' 


_p_ 

Poo 


p m 



T’ = 



(2.14) 


t 

1 = 0 V 2 

roc * og 

where the nondimensional variables are denoted with an asterisk, t m is the nondimen- 
sionalized time, and x‘ are the nondimensionalized coordinates. Applying this nondimen- 
sionalization to the compressible Navier-Stokes equations, one obtains a set of equations 
which assume essentially the same form as the original formulas provided that new material 
constants are defined as follows: 


P 


j£_ 
Re i 



(2.15) 


c v = b(7-i )Mir 1 

In equation (2.15), Moo is the free stream Mach number which is given by the relation 


Moo = , 00 - (2-16) 

V 7(7 - IKToo 

Throughout this report, asterisks and tildes are omitted and it is understood that the nondi- 
mensionalized forms of the governing equations are in use. 

This section is concluded with a list of some nondimensional quantities of interest: 
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• pressure coefficient 


c P = 


P- Poo 


0.5 p 


1/2 

CO * oo 


= 2(p- - 




(2.17) 


• skin friction coefficient 


C J 



(2.18) 


where the nondimensional shear stress r* is calculated using the nondimensional form 
of (2.8) (with p and A replacing p and A, respectively) 


• heat flux coefficient 


where 


Ch = 


0.5 pocV*, 


= 2q' 


(2.19) 


q = q t ni and q’ = q“ n, (2.20) 

and where n = (n,) denotes a unit normal vector. The nondimensional heat flux 
vector q’ is evaluated using formula (2.9), the nondimensional temperature T* and the 
coefficient of conductivity 


_ ICvP 
Pr 


( 2 . 21 ) 


2.2 Taylor-Galerkin Algorithms 

This section presents two Taylor-Galerkin algorithms for integrating the Navier-Stokes equa- 
tions forward in time. The first algorithm is a implicit/explicit one-step approach where all 
of the components of the solution are handled simultaneously. The second approach is a 
rather unique two-step approach whereby the Navier-Stokes equations have been split into a 
convective part and a diffusive part and are subsequently combined to advance the solution 
in time. 
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2.2.1 One-Step Taylor-Galerkin Algorithm 


A general family of implicit Taylor-Galerkin methods is summarized in this section. These 
methods can be made explicit or implicit by the appropriate selection of implicit parameters 
(see below). The details concerning the theoretical formulation can be found in [27] and 
references therein. 

Given a domain f l C the compressible Navier-Stokes equations are characterized by 
a system of conservation laws of the form 


u, t + Ff ti = Fl xen,t> o 

(2.22) 

accompanied by an initial condition 


17(*,0) = U 0 (x) 16(1 

(2.23) 

and appropriate boundary conditions. 



Second-Order Taylor Expansion in Time 

As a starting point, let us assume that the solution U n at time step t n is known and the 
solution £/ n+1 at time t n+l is to be calculated. Formally, the values of the solution at times 
t n and time t n+1 can be expressed by the second-order Taylor series expansion about an 
arbitrary time t n+a where a is an implicitness parameter with values between zero and one: 


At 2 

t/ n+1 = U n+a + (1 - a)AtU™ +a + (1 - a) 2 — U n u +a -I- 0(A< 3 ) 

At 2 

U n = U n+a - aAtU"+° + a 2 =—U^ t + ° -I- 0(At 3 ) 

£ 


(2.24) 


By subtracting these two formulas a formula is obtained for an increment of the solution 
between steps n and n + 1 : 


AU = U n+1 - U n = AtU” + ° + (1 - 2a) ~U n t ^ Q + 0(At 3 ) (2.25) 

Observing that: 

L/” +a = L7” +/3 + 0((a - 0)At) (2.26) 

a second implicitness parameter can be introduced into equation (2.25) while still preserving 
the second-order accuracy, 

A j2 

AU = AtU" + ° + (1 - 2q)—U^ 0 + 0(At 3 ) 

£ 


(2.27) 



The next step is to express the quantities evaluated at times < n+ ° and t n+0 by quantities 
evaluated at the basic steps t n and < n+1 , 

C/” +a = l/" + aAU t + 0(At 2 ) (2.28) 

Uu +0 = U* + 0Al/« 4- 0(At 2 ) (2.29) 

Substituting these formulas into equation (2.27) yields a two-parameter expansion: 

A / 2 

At/ = At {U n t + oAt/j) 4(1- 2a)— (t/" + /?At/«) + 0{At 3 ) (2.30) 

Now, following the procedure outlined by Lax and WendrofF [21], the governing equations 
can be used to replace time derivatives with spatial derivatives. This substitution yields a 
formula for the first derivatives: 


U, = F l - F 




and for the second-order derivatives: 


l/« = 


where Ai, Pi , and it,, are the Jacobian fluxes 


(2.31) 


R„ {F' u - Fi t ) . + P, (pr,* - Ff,.)] . - [* {rL - *£.)] . (2.32) 


A,= 


dF * 


Pi = 


dF 


V 


R>: = 


dF) 


(2.33) 


dU “ ' dU dUj 

The convective terms in equation (2.32) involve spatial derivatives up to fourth order. Lim- 
iting this formula to terms with second-order derivatives, which can be effectively handled 
by C° continuous finite elements, yields the following approximation for the second-order 
time derivative 


V„ = (AiFfj,) . + 0(/i, it) (2.34) 

where 0(n,k) represents a quantity of the order of the viscosity parameters in the Navier- 
Stokes equations. Substitution of formulas (2.31) and (2.34) into the incremental equations 
(2.30) gives the implicit formula: 

AU = At - Ff;) + Q (AF£ - AF£)] 

+ (1 - 2a)^ [(AFy . + 0A (AFy ,.] + o (A! 3 ) + 0(n,k)0 (AO) ( 
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For the sake of maximum generality, a third implicitness parameter can be introduced 
by observing that 

aAFj; - aAF? = 7 A F?- - qA Ff- + (a - 7 )0(/x, k)O(At) (2.36) 

Substituting this expression into (2.35) yields the three-parameter implicit form for the 
increments of the conservation vector 


A U = to l(FV* - Ff;) + -rAFji - orAFg] 

(2.37) 

+ (1 - 2a)^- [(^iFf fc + 0A (AiFtJ .)] 

Equation (2.37) represents a nonlinear formula for increments of the solution U at a 
given time step. This formula is nonlinear due to nonlinear dependence of the fluxes and 
Jacobians on the solution U . This equation can be linearized, while still preserving second- 
order accuracy, with the resulting incremental formula: 


AU + a At (A"AU) t - -,At [(KJAt/.,) . + (P.'AI/) ,] 

- (1 - 2c)l3 (A'A'AU,) . (2.38) 

= to ( Fj; - Ff;) + (i - 2a) ^ 

The particular form of this evolution equation used in the finite element approximation 
corresponds to setting a = 0 in (2.38). For this special case one obtains 


AU 


7 Ai [(«" AU,) + (P^AUl 

' ' 1* , 



(2.39) 


Weak Formulation 

In order to obtain a weak variational formulation of the incremental equation (2.39), we 
introduce the space of test functions 


13 



W = {v = (Vi, V 2 . . . V M ) s.t. V, e H\n) and K = 0 on T c } (2.40) 

where M is the number of conservation variables, H*(£l) is the usual Sobolev space of func- 
tions with derivatives in T 2 (fi), and Tp is the boundary with specified Dirichlet boundary 
conditions. After multiplication of the incremental equation (2.39) by an arbitrary test func- 
tion V^x) € W 7 , integrating over the domain ft and application of the divergence theorem, 
the following weak formulation of the problem is obtained: 

Find AU € W s.t. V V <E W: 



AUV + 7 A tR^AUj ■ V,,- 


+ 7 AtP?AU-Vj + 



A*A]AU' } ■ V 


0 


dn 


-[ ( 7 A in t RlAU <r V 

Jan \ 

At 2 \ 

+ 7 A tn,P?AU ■ V + fi — mA^AJAUj • VjdS 

= J n ( A( (J!f ■ - /-!'") • V,, - a:a;u, ■ v. t j 

+ / 8n (-Ain, (Ff" - FV») • V + ^ • v) <iS 


(2.41) 


It can be shown by the selection of appropriate test functions that the solution of this 
problem is, in the sense of distributions, the solution of the boundary-value problem, to- 
gether with appropriate boundary conditions. Additional details concerning implicit Taylor- 
Galerkin methods may be found in references [11,15,22,31]. 


2.2.2 The Two-Step Algorithm 

A second algorithm investigated during the course of this project for solving Navier-Stokes 
equations is based on a two-step approach [8,10,28]. The method consists of advancing 
the solution in time by performing interchangably two steps associated with the convection 
operator E and the diffusion operator H corresponding to inviscid and viscous terms in 
equation (2.1): 


I/ n+1 = G{t)U n (2.42) 

where G(t) = H(t)E{t). The convection operator E{ t) is defined by: 
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(2.43) 


(E(t)U 0 )(*)~U(x,i) 
where U(x,t) is a solution to Euler equation: 

’ V, + Ff (U)i = 0 

i 

. U(x,0) = U o (x). 

The diffusion operator H(t ) is defined by: 

(JT(t)t/ 0 )(x) = f U(x,t) 

where U(x,t) is a solution to: 


(2.44) 


(2.45) 


' U, t = /?'(£/),■ 

. U(x,0) = U o (x). 


(2.46) 


Problems (2.44) and (2.46) must be augmented by appropriate boundary conditions, a 
detailed discussion of which is given in the following section. It should be mentioned that a 
different composition of operators H(t) and E(t) gives a three-step procedure of the form 


G(t) = H(t)B(l)H(|) (2.47) 

which is second order accurate while our two-step procedure is only first order accurate 
in time. This however is not of our concern since we are interested only in steady-state 
solutions. 

In the numerical implementations, exact operators E(t) and H(t) are replaced by their 
discrete approximations. The motivation of applying this two-step approach to solving 
Navier-Stokes equation was that different solvers or even different spatial approximations 
could be used for the Euler and viscous steps. We may take advantage of this flexibil- 
ity especially in modeling boundary phenomena: solving, for instance, the convection step 
with a specialized and very efficient Euler solver and linear approximation, while using very 
accurate higher order approximation in boundary layer zones in the viscous step. 


The Euler (Convection) Step 

In this work, we approximate the Euler step (2.44) by a second order Taylor- Galerkin scheme 
while the viscous step (2.46) is approximated by a first order finite difference approximation. 
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This reduces the problem to solving at each time step elliptic-like boundary value problems 
for which we employ our h-p finite element method for the spatial approximation. 

A Taylor- Galerkin scheme for the Euler equations is obtained by simply setting F]' and 
the various derivatives of F J to zero in equations (2.39) and (2.41). For completeness, this 
procedure is summarized below. 

Given a domain Q C , the compressible Euler equations are characterized by a system 
of conservation laws of the form 


U it + Ff (U),i = 0 a: € fi,t > 0 (2.48) 

accompanied by an initial condition 

U(x,0) = U o (x) xzSl (2.49) 

and by appropriate boundary conditions. 

Now taking equation (2.39) and setting R,j, P, and F \ to zero one obtains an evolution 
equation for the conservation variables 


AU - 0 ~ (AWAUj) . = 

- AtFf; + ^ (4”Fg) . 

The corresponding weak formulation follows from (2.41) as 
Find AU € W s.t. VV G W 

J a (au ■ V +0^-a:a;auj ■ v.^j d(i 
-J K (0^^-v)0s 

= LtF?" ■ V, - . V,) in 

+ J sn f-Aln.Ff" ■ V + $£n i A?A?Ujv') iS 


(2.50) 


(2.51) 


The Taylor-Galerkin method for solving the Euler step is unconditionally stable for 
ft > 0.5 independently of the approximation in the space variables. The second-order terms 
present on the left-hand side modify the L 2 -projection and contribute to the stabilization of 
the method. Additional details are given in [8]. 
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The Viscous (Diffusion) Step 


The second step in the two-step method is the viscous or diffusion step defined by equation 
(2.46). In this step, the density remains unchanged and evolution equations for the mo- 
mentum and energy may be fully decoupled provided that the boundary condition for the 
momentum equations can be formulated so that they do not contain energy terms. Under 
this assumption, one arrives at a system of two equations to be solved simultaneously for 
the momentum components 


drrii 

~dt ~ Ti} ' j 

and a single scalar equation to solve for the total energy 


(2.52) 


de 

~di 


( T uU,)„ + 


(2.53) 


As a starting point in the solution of the momentum and energy equations, we introduce 
a Taylor series expansion of the conservation vector U about an arbitrary time t + 0At 


U{t + At) - 0AtU ti (i + At) = U{t) + (1 - 0)AtU, t (t) + 0{At 2 ) (2.54) 

where j3 is a constant between zero and one. Again using the Lax-Wendroff procedure for 
replacing the time derivatives by spatial derivatives one obtains 


and 


m] +l - 0AtT”+ l = m" + (1 - 0)Atr^ 


(2.55) 


e "+i _ 0A t £ (t" +, u" +1 + kT" +1 ) . 

= e" + (l-m(E(r;u" + <c7l) i 

.= 1 


(2.56) 


The procedure for solving the system of equations then becomes: solve the inviscid step 
using the implicit Taylor-Galerkin method outlined in (2.51), next solve the momentum 
step defined by equation (2.55), and finally solve the energy equation for the temperature. 
Combining these results as indicated by equation (2.42) advances the solution of the Navier- 
Stokes equation found in time by At. 
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The variational formulation for the momentum and energy equations are obtained using 
exactly the same procedure that was used for the convective or inviscid step. Performing 
these steps one arrives at the following variational formulations: 

Find m” +1 such that 


/ m’ l+1 v'<m + /?A< ( ^ + 1 v itj dn 

Jn J n 

- f3At [ T^njVidS = / m^VidQ 

Jd n } Jn 

- (1 - /3)At [ T?jVijdn + (1 - p)At [ T^njVi 

j n j ovi 


(2.57) 


for every K 


and 


Find e" +1 such that 


J e n+1 VdU + /3 At jf KT? +l V,dn 

~/3At I KT" + 'n t Vds = 

Jan 

J e n Vi(l - (1 - P)At jf kT"V,iKI 

+ (l-/?)Af J KT^n.ViS 

+ At rj+'uj*' + (1 - ?)r^V,dn 


(2.58) 


for every V 

Note that in the variational formulation of the energy equation one obtains a volume 
integral which is a function of the stress and velocity at time t + At. By allowing the 
viscosity parameters A and fx to lag a time step, we are able to first solve the momentum 
equation for the velocity components and then explicitly evaluate this final term. 


2.3 Boundary Conditions 

This section presents a general overview of COMCO’s approach to prescribing boundary 
conditions for the Navier-Stokes equations. This approach is based upon a linearized stabil- 
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where ii' = p 


XT + V 


,2 


Figure 2.1: The symmetrizer. 

ity analysis of the Navier-Stokes equations which results in the following entropy stability 
condition that must be satisfied on the boundary 

f ( -6U T A 0 A n 6U - SU T A 0 St n ) ds > 0 (2.59) 

Jan \2 J 

In this equation, SU is the variation in the conservation vector, A n is the normal Jacobian 
matrix defined by A„ = A,n,, A 0 is the symmetrizer of the Navier-Stokes equations shown 
in Fig. 2.1, and St„ is the variation in the normal viscous flux. Additional information on 
the symmetrizer and stability conditions (2.61) can be found in references [8,9]. 

2.3.1 Boundary Condition for the One-Step Algorithm 

We begin our discussion of boundary conditions for the implicit one-step Taylor- Galerkin 
methods by quoting a result of Strikverda [29], which specifies the number of boundary 
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conditions necessary for well-posedness of the linearized Euler and Navier- Stokes equations. 
These results are summarized for two-dimensional problems in Table 2.1. 

Table 2.1 


Type of 
Boundary 

Euler 

Navier-Stokes 

supersonic inflow 

4 ess 

4 ess 

subsonic inflow 

3 ess 

3 ess + 1 nat 

subsonic outflow 

1 ess 

1 ess -1- 2 nat 

supersonic outflow 

0 

0 ess -f 3 nat 

no-flow 

1 ess 

1 ess -1- 2 nat 

solid wall 



— isothermal 

— 

3 ess 

— heat flux 


2 ess -f 1 nat 


In this table, “ess” denotes the essential boundary condition and “nat” denotes natu- 
ral boundary conditions. The essential conditions are to be imposed on the characteristic 
variables rather than the conservation variables. It is important to note that the numbers 
presented in the table are true for problems that are not regularized. If artificial diffusion is 
built into the algorithm or added explicitly, natural boundary conditions should be imposed 
on these terms even for Euler problems. Moreover, since artificial diffusion can affect all 
conservation variables, the number of natural boundary conditions for these terms should 
actually be one more than for the (nonregularized) Navier-Stokes equations. 

Before launching into a full discussion of the various classes of boundary conditions, it is 
useful to first recast the boundary integrals in terms of the characteristic variables (Riemann 
invariants). In the variational formulation a typical boundary term is of the form 

AimAU • V = A n AU • V (2.60) 

where A n is a nonsvmmetric matrix. The matrix A n can formally be represented with 
respect to its own eigenbasis as 

M 

A„ = Y. A .(c.®K) ( 2 - 61 ) 

cr=l 
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where b Q and c a are the left and right eigenvectors, respectively. The eigenvalues A a for the 
two-dimensional case are 

A, = u n c 

A 2 = Un 

A3 = U„ 

A 4 = U„ + c 

where u n is the velocity normal to the boundary and c is the speed of sound. The expressions 
for eigenvectors b Q and c a can be found in references [8,12,29]. Note that a positive value of 
A 0 means that the corresponding characteristic exits the domain across a given boundary, 
while negative values of A 0 correspond to signals entering the domain. As a general rule, the 
characteristic variables corresponding to characteristics entering the domain (negative A Q ) 
need to be specified as the essential boundary conditions, while the characteristic variables 
exiting the domain are continued across the boundary from the interior. 

The characteristic variables A U 0 are defined as components of the vector AU in the 
eigenbasis of A n so that 


(2.62) 


AU = (AU ■ b 0 )c Q = AU a c 0 
V = (V-c a )b a = V a b a 

With these definitions, the boundary formula (2.60) can be presented in terms of character- 
istic variables as: 

M 

A n AUV = Y,KV a AU a (2.63) 

0=1 

The above representation is very useful in the formulation of essential boundary conditions. 

In the following sections the formulation of the boundary conditions for various boundary 
types is discussed. Additional details can be found in [12,13,25] and references therein. 


Supersonic Inflow 


On a supersonic inflow boundary, the values of all the characteristic variables (thus also of 
all the conservation variables) are specified as the upstream values. Formally, this means 
that 


U = U 


(2.64) 


or, in incremental form, 


AU = U - U n on dQj 
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In practical applications these conditions are enforced by the penalty method, which is 
obtained by adding to the variational formulation a term 


- [ \ AU -(U- l/ n )l • Vds 

t JdUj 1 v /J 


e being a small parameter. 


Supersonic Outflow 

On a supersonic outflow boundary, all characteristic variables propagate along characteristics 
so as to be continued from the interior of the domain and no essential boundary conditions 
are imposed. The explicit algorithm application of supersonic outflow boundary conditions is 
achieved simply by calculation of the boundary integrals. In the implicit algorithm, natural 
boundary conditions are imposed on viscous and second-order terms. 

Natural boundary conditions on viscous terms are imposed by observing that the viscous 
boundary terms on the left-hand side of the variational equation (2.41) can be interpreted 
as 

K)'jAUj = [ + PiTuViVj) AUjds 

J dQ. 0 


= / -7 A iAF^rfs 

Jd Q 0 

where the components of AF„ are 

AF\ = {0, Aa ln , A<7 2n , Aq n } 7 

and )P(x) are the shape functions. In order to formally impose natural boundary conditions, 
the above terms are transferred to the right-hand side with prescribed values of AF^, so 
that the new right-hand side is 

= fAtAFl^ids 

JdUo 

Note that since the mass flux due to viscous terms is identically zero, this procedure actually 
imposes only three natural boundary conditions (in two dimensions). 

The choice of the actual conditions is somewhat arbitrary. Currently two options are 
implemented, namely, 

• zero change of flux at the time step (frozen viscous flux): 

3F^ = 0 
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• zero total viscous flux, enforced by: 


AF V n = 0-F v n n 

Obviously, on an outflow boundary adjacent to a solid wall the viscous fluxes are not zero 
and the first procedure is more appropriate. 

In order to ensure well-posedness of the problem, proper natural boundary conditions 
should also be imposed on the second-order terms. Analogously as in the viscous case, 
second-order terms on the left-hand side of equation (2.41) can be transformed to the form: 

t A 

KfjAUj = 1-0 -=-*/ (AnAFjj) ds 

J dQo 

This term has no simple physical interpretation and therefore the selection of the natural 
boundary condition to be imposed is somewhat difficult. The procedure adopted by Demkow- 
icz, Oden, and Rachowicz [8] is to decompose the above term into components normal and 
tangential to the boundary and impose boundary conditions only on the normal term (sym- 
metry boundary conditions). In this work, a slightly different procedure is applied, according 
to which the whole term is transferred to the right-hand side with certain prescribed values. 
This corresponds to imposing natural boundary conditions on 

A n AFjj = A n AjAUj 

The actual boundary condition applied is to set the total value of this term to zero, so that 
the prescribed value at the time step is 

AnAjAUj = 0-A n AjUl (2.65) 

An interpretation of this condition can be obtained by observing that for Euler problems, 
the enforced condition is A n U = 0, or in terms of characteristic variables, 

^ a i r Q b a = 0 o = l,...,M 

Since, on the supersonic outflow, all the eigenvalues satisfy A 0 > 0, this condition means 
that the characteristic variables do not change in time as they exit the interior across the 
boundary. 

A somewhat more appealing interpretation can be presented for the simple two-dimen- 
sional advection equation 

U = a, l/,, * = 1,2 

for which the characteristics are straight lines defined in space-time by the vector c = 
{ai,a 2 ,l}. The natural boundary condition corresponding to (2.65) is: 

nididjUj = 0 
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or equivalently: 


a n (DU , a) = 0 

where (DU, a) denotes the directional derivative of u in the direction of a. This means 
that on the outflow boundary ( a n > 0) U must be constant along advection lines. 

On a contour plot, this forces contours of U to be parallel to advection lines on the bound- 
ary. It can be observed that the condition applied by Demkowicz, Oden, and Rachowicz [8] 
causes derivatives of U along the normal to the boundary n to be zero, which corresponds 
to contour lines normal to the boundary. 

Subsonic Inflow and Outflow 

On a subsonic inflow or outflow’ boundary, essential boundary conditions are imposed only 
on the characteristic variables corresponding to negative eigenvalues X a . For each of these 
terms, the corresponding condition is 

A U a = AUb a 

= (U-U n )-b Q 

where the prescribed far-field values of the conservation variables are denoted by U. The 
penalty term enforcing essential boundary conditions on the increments of selected charac- 
teristic variables A U 0 is of the form 

Kjj = -/ (c a ®b 0 )V } Vjds 
e Je n 

Rj = - I \(U - U n ) • 6 J c a tf/<fc 
t JdQ L J 

Nodewise application of these conditions is obtained by replacing the shape functions with 
Dirac delta functions associated with boundary nodes. 

For the characteristic variables with nonnegative eigenvalues, continuation from interior 
conditions are employed. These conditions for selected characteristic variables involve rather 
complicated formulas. Thus, for practical applications it is better to observe that, since 
the penalty procedure actually overrides any other conditions, the continuation condition 
can first be applied to the whole vector of conservation variables (by the supersonic outflow 
procedure), and then the above penalty method can be applied to selectively enforce essential 
boundary conditions. 

In practical implementations of subsonic outflow boundary conditions, some authors 
choose to prescribe a value of pressure p rather than the value of the characteristic vari- 
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able. The incremental form of the pressure boundary condition is 


Ap = p-p n 

and can be enforced by a penalty method. The corresponding penalty term in the variational 
formulation is 

i[Ap-(p-p”)]-<(Ap) (2.66) 

The pressure increment A p can be expressed in terms of conservation variables as 

Ap = % AU = dAU 

T 

where d = (7 — 1){^S Ui,u 2 ,lj (in the two-dimensional case). Therefore, the penalty 
term becomes 

l -\d-AU-{p-p n )){d- V) (2.67) 

where V is a test function. The corresponding stiffness matrix and right-hand side are then 

K,j = - / d®dVjyjd$ 

R, = - [ d(p-p n )' i it]ds 

e Jd n 

Solid Wall 

There are basically two types of solid wall boundaries: 

• adiabatic walls with prescribed zero heat flux (M-th component of the viscous flux 
vector): 

= ^n(A/) = 0 

• isothermal walls with specified temperature: 

T = T 


In addition to the above conditions, zero velocity (zero momentum) conditions are also 
specified on the solid wall. These conditions are easily enforced by the selective application 
of a penalty method, similar to the supersonic inflow procedure. The incremental form of 
the adiabatic condition of zero heat flux is 

4JV, = -F (2-68) 
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This natural boundary condition is applied by formally transferring the viscous terms cor- 
responding to the energy equation from the left-hand side of the variational equation (2.41) 
to the right-hand side and setting the increment of the heat flux according to (2.68). It is 
of interest to note that since the viscous terms do not directly affect mass fluxes and the 
momentum equations are overridden by the penalty method, in practice all viscous contri- 
butions can be skipped on the left-hand side. 

On the isothermal wall the additional boundary condition is a prescribed temperature T 
or, equivalently, a prescribed specific energy e. Since the kinetic energy is zero on the wall, 
the above condition can be expressed in terms of conservation variables as: 

e 

- = e 

P 

In the incremental form this condition becomes 

-A e - = (e - e") (2.69) 

PP 7 

This condition is imposed via the penalty method. It should be noted that there are available 
a variety of possible forms of the penalty terms, depending on the form of the test term 
applied to condition (2.69). One possibility, which appears to be the most natural and yields 
a symmetric contribution to the stiffness matrix, is obtained by testing equation (2.69) with 
its own variation: 

7 [(^ Ae ~ J^ p ) - (l ~ £n) ] ' G V,M > ” 

where Vjq denotes a selected component of a test vector: V(i) for density and I \m) fo f energy. 

This approach leads to two penalty conditions affecting both the continuity and energy 
equations. Therefore it is not in agreement with the physical situation because, while the 
solid wall can supply heat to maintain a prescribed termperature, it cannot supply mass 
for this purpose. For this reason, another form of the penalty term should be used, which 
enforces a prescribed temperature by altering the energy equation only: 

i f ( -Ae - -% aA - (e - e n ) ■ V,,,, 

* \P P J J 

The corresponding terms in the stiffness matrix and the right-hand side are 


Kjj = / kVjWjds 

J W 

Rj = / r'P/ds 

J dQw 


(2.70) 
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where the kernel matrices k and r are defined (in two dimensions) as: 


k 


0 



0 0 0 ' 
0 0 0 
0 0 0 
0 0 1 


r 


1 

e 


P 

0 

0 

0 


p{e - e n ) J 


Again, nodewise enforcement of these conditions can be obtained by replacing shape func- 
tions with Dirac delta functions. 


For the regularized problem some additional artificial terms (fluxes) occur on the solid 
wall due to second-order terms and explicit artificial dissipation. These fluxes are forced 
to be zero by means of natural boundary conditions, in the same manner as the supersonic 
outflow. 


No-Flow 

The basic condition of the no-flow boundary is that the normal velocity is zero or, equiva- 
lently, that the normal momentum is zero, 

m,n, = 0 


In the incremental form this becomes 


Am, 7i, = — m"n,- 

This condition is easily enforced by a penalty function, with the addition of the term 

-/ (Am.nj + m "n,-) V {l+i) ds (2.71) 

e Jau w 

in the variational formulation, where V(i+,) is the component of a test function corresponding 
to momentum m, . The resulting stiffness matrices and right-hand sides are of the standard 
form (2.70), with kernels (in two dimensions): 

0 0 0 0 ' 

0 ni«2 0 

0 ri2ni 02^2 0 

0 0 0 0 
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r = 

For viscous flow, the additional conditions on the no-flow boundary are the two natural 
boundary conditions: 

o n > = 0 

9n = 0 

where a ni is the skin friction on the boundary and q n is the normal heat flux. The appli- 
cation of these natural boundary conditions follows the procedure discussed in preceding 
subsections. Similarly, as on the solid wall, all artificial fluxes are forced to be zero on the 
no-flow boundary. 

2.3.2 Boundary Conditions for the Two-Step Method 

This section presents a brief overview of the boundary conditions for the two— step algorithm 
outlined in Section 2.2.2. In general, the boundary condition can be constructed separately 
for the convection (Euler) step and the diffusion (viscous) step. This holds provided that 
they guarantee stability of those steps and possess appropriate asymptotic properties as 
the viscosity constants approach zero ( Rei — ► oo). With this in mind, we summarize the 
different classes of boundary conditions appropriate for the Euler step, momentum step, and 
energy step. 


0 
n i 
n 2 
0 


Euler Step 


The boundary conditions for the Euler step in general follow directly from the one-step 
algorithm with the following special conditions. 


1. Contributions of viscous fluxes to boundary terms in variational equation (2.41) are 
omitted. 

2. The essential boundary condition on a solid wall is limited to enforcing a zero normal 
component of the momentum vector. It is accomplished by means of a penalty method, 
i.e., by adding the following contributions to the stiffness matrix: 


1 

£ 


f (JJ2 n x ^3 n y) {&U2 n x "h &U3 n y) dS 

JdQ 


where e is a small penalty parameter. 


28 



Boundary Conditions for the Momentum Step 

Boundary conditions for the momentum and energy steps were constructed such that the 
viscous term in expression (2.59) results in a positive production of entropy and the resulting 
boundary terms in boundary value problem (2.57) and (2.58) make these problems well posed. 
These boundary conditions can be listed as follows: 


Case 1 


Case 2 


Case 3 


Case 4 


Open Boundary — Supersonic Inflow 

Full Dirichlet boundary conditions are prescribed, 


m" +1 = m’ n 


(2.72) 


where mj n are the momentum components of the same inflow vector as used in 
the supersonic inflow boundary conditions for the Euler step. 

Open boundary — Subsonic Inflow 

The same Dirichlet boundary conditions are used, but with the inflow vector 
replaced with the solution from the Euler step, i.e., 


. n+1 

m. = m. 


Open Boundary — Subsonic Outflow 

Mixed boundary conditions are used, 


m ” +1 = m n 


• 1 _ y ^ 

*3 3 


where m n is the normal component of the momentum, 

m n = mjTij + m 2 n 2 

and t, is the tangential viscous stress vector component, 

t„ = (r 22 - Ti^nw + r 12 (nj - n 2 ) 

Open Boundary — Supersonic Outflow 

Full Neuman boundary conditions are applied, 


(2.73) 


(2.74) 


(2.75) 


(2.76) 


T n ."t^n ■ — T— 71 • 
T ij 11 1 'tjJ 


(2.77) 
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Case 5 


Solid Wall Boundary Conditions 

Full Dirichlet boundary conditions are used, 

mr +1 = 0 (2.78) 


Case 6 Symmetry Boundary Conditions (of the first kind) 

Mixed boundary conditions are applied, 


where m n 
vector. 


m 


n+l _ 0 


dm ?* 1 

dn 


(2.79) 


= 0 


and m, are the normal and tangential components of the momentum 


Case 7 Symmetry Boundary Conditions (of the second kind) 

Full Neuman boundary conditions are applied: 

dm ?-' = o (2.80) 

dn 

As in the case of Euler equations, substituting (2.80) into the boundary integral 
in (2.57), results in some non-zero terms which must be included in the stiffness 
matrix calculations. 


In the finite element code, all of the essential boundary conditions have been implemented 
using the penalty approach, i.e., replacing the full Dirichlet boundary conditions with 

+ e0At r," +1 nj = m‘" (2.81) 

j=i 

and the first of (2.74) conditions with 

m n n +1 + epA tr; +1 = m” (2.82) 

where r" +1 is the normal viscous stress 

T n = 2 T 0 n ' n J ^ 2 - 83 ^ 

i,j=l 
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Boundary Conditions for the Energy Step 


Case 1 Temperature (energy) Prescribed 

A single Dirichlet boundary condition is applied 


e 


n+l 


e 


( 2 . 84 ) 


The choice of e will vary with the particular kind of boundary: 


• for supersonic inflow e corresponds to the inflow vector. 

• for subsonic inflow e is evaluated using the solution from the previous step. 

• for a solid wall with temperature prescribed e corresponds to the prescribed 
temperature on the wall and the density p (unchanged in the viscous step). 

Note that the two-step method eliminates the contradictions resulting from the 
discussion of the solid wall boundary conditions with temperature prescribed, as 
the density remains unchanged (see [10]). 


Case 2 


Heat Flux Prescribed 

A single Neumann boundary condition is applied: 


dT 


•n + l 


dn 


= 9 


The heat flux q is calculated in the following way: 


( 2 . 85 ) 


• for subsonic and supersonic outflow, q is evaluated using the solution from 
the previous step 

• for an adiabatic wall and symmetry boundary conditions of both kinds, q is 
assumed to be zero. 


2.4 Artificial Viscosity 


In order to suppress spurious oscillation of the solution an artificial dissipation is introduced 
as an additional flux in the Navier-Stokes equations in the form 


Ut + FZ-Fl + Fi, ( 2 . 86 ) 

where denotes the artificial dissipation flux with corresponding Jacobians defined as 

dFf 


P? = 


dU 


r a _ dFf 

3 au 0 


31 



The advantage of this approach is that the artificial dissipation can be treated using exactly 
the same formulation and procedures as for the natural viscosity. In the one-step algorithm, 
for the sake of generality, a fourth implicitness parameter 7 is introduced for the terms 
associated with the artificial dissipation. In the calculation of the stiffness matrices, right- 
hand sides and boundary terms, the same formulas are used as for the natural viscosity. 
Similarly, for the two-step algorithm, the artificial viscosity is implemented as an additional 
“viscous” flux when solving the equation of the Euler equation. 

In this work, two commonly used forms of artificial viscosity have been studied and 
implemented as follows: 


• the classical Lapidus viscosity [20] 

F A = kuU,i (2.87) 

with 

k i{ = c k h 2 |u t ,,| (2-88) 

The Jacobians P A and R A can be defined by a straightforward differentiation of (2.86). 


• the generalized Lapidus viscosity due to Lohner, et al. [23] 


F a = Ik 


dU 

dt 


or F A = kUl } Uj 


(2.89) 


with 


k = Ckh 2 (£ ■ V(u • £)) 
£ = VM/|V|u|| 


(2.90) 


where h is the element size, c k is an arbitrary coefficient controlling the amount of 
dissipation, k is a solution dependent set of coefficients, / is a unit vector parallel to 
V[u|, and u is the velocity vector. The tensor product IJj ensures that the artificial 
viscosity acts in the direction normal to the shocks. The Jacobians P A and R can 
be defined by differentiation of (2.88). For simplicity, dependence of k and t on the 
solution is disregarded, in the definition of jacobians so that 



R* = kUl,I , 


where I is the identity matrix. 


The drawback of the above formulations is that for elements with high aspect ratios the 
element size h is not a clearly defined quantity. The wrong choice of h may cause oscillations 
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if h is too small or considerable smearing of shocks if h is too large. We have developed 
a modification of the generalized Lapidus viscosity (2.88) such that it uses more precise 
information about the geometry of the element than just h but still preserves its original 
form for square elements. In terms of the contribution to the element stiffness matrix, the 
generalized Lapidus viscosity can be written in a slightly different form as 


c k At 



v,, k i} Uj<m 


where fc.j = and k is a solution dependent scalar 


k = 


d(u ■ t ) 

dt 


(2.91) 


(2.92) 


The basic idea of constructing the modified artificial vsicosity is to perform the calculation 
of (2.90) on the master element, which has a fixed size ( h = 2), then map it back to the 
physical element domain. As a result of this procedure the modified artificial viscosity is of 
the form 


C/t At h 2 



(2.93) 


where k,j is given as 


- _ chj (hj 

0 ' at, kpq at 


(2.94) 


»P “M 

and k pq is of the same form as in (2.89) except that the unit vector l is taken as 


< = V ( |r|/ |V<M| 


(2.95) 


to make the viscosity work in the direction perpendicular to the shocks on the master element. 
V f denotes the gradient calculated on the master element coordinate. It can be shown that 
for square element this modified artificial viscosity does coinside with the original expression 
(2.90). We have applied this modified viscosity in several test problems where element 
with high aspect ratios are used, and have found it to be more effective than the artificial 
viscosities defined by equations (2.88) and (2.90). 
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3 An h-p Finite Element Method 

The aim of adaptive methods in computational fluid dynamics is to optimize the compu- 
tational process: to obtain the best results for the least effort. The cost functional in this 
optimization process is the numerical error measured in some appropriate norm, both the 
global error over the entire computational domain and the local error over each gridcell. 
The central parameters are the conventional mesh parameters that govern local accuracy: 
the mesh size h, the order p (e.g., the spectral order) of the local approximation, and the 
location of gridpoints. 

In the h-p FEM one can control both the local mesh size and spectral order of approxi- 
mation simultaneously. Such a flexibility allows one to distribute degrees of freedom in an 
optimal way: a large density of degrees of freedom can be used in computational regions 
with very irregular behavior of the solution while a relatively rough approximation is used 
in subdomains where the solution is smooth. This suggests that the h-p method may use 
an optimal number of degrees of freedom to achieve a prescribed accuracy. In addition, 
recent work in the area of approximation theory [2,3] suggests that an extra gain in accuracy 
can be obtained if the enrichment of the mesh is performed in two combined ways: first by 
reducing the size of elements h and second by increasing their spectral order p. The problem 
of how to combine these two kinds of refinements so that the improvement in accuracy is 
the best possible is a very complex issue. In general, its strict mathematical solution is not 
known, however, there exists heuristical knowledge on the use of h-p FEMs for many classes 
of problems. 

In practice the reduction of the mesh size h can be achieved in two ways: by subdividing 
elements into smaller sons, or by so-called remeshing, i.e., generating a completely new 
mesh with a given distribution of h. Our implementation of the h-p FE method uses the 
first approach. We break two-dimensional quadrilateral elements into 4 element sons and 
three-dimensional hexagonal elements into 8 sons. 

The success of such a complex adaptive scheme depends upon several properties of the 
adaptive process: the data structure, the adaptive strategy, the techniques for a posteriori 
error estimation, and the flow solver. The potential payoff of a successful h-p adaptive strat- 
egy is substantial: exponential rates of convergence may be attained, meaning that complex 
flow features can be resolved using orders-of-magnitude fewer unknowns than required by 
conventional methods. 
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3.1 A Variational Formulation 


Let us first specify the class of problems to be solved with an h-p FEM. In this development, 
we follow the detailed presentation in our papers [9, 11, 25]. 

Let Cl be an open bounded domain in R n , n = 2,3 with a sufficiently regular boundary 
dCl. In what follows, we shall restrict ourselves to a class of problems that can be formulated 
in the following abstract form: 



Find u € X such that 1 

B(u,v) = L(v) Vv € X J 

(3.1) 

Here: 

X = X x X ■ • • x X(m times) 

(3.2) 

where X a subspace of the Sobolev space of first order, B( •,•) is 

X x X of the following form 

a bilinear form on 


m 

B{u,v) = ^ Bij(ui,Vj) 

i t j=l 

(3.3) 

where B tJ ( • 

, •) are bilinear forms of scalar-valued arguments of the type 



B,j(u,v) 



/ ( JL* . , du dv ^ k du 1 

= h^^ + £-’^ v+c - ,uv \ 

(3.4) 


+ / d,,uv ds 
JdU 


and 

L(-) is a linear form on X 

(3.5) 

of the form 

L{v) = Y.Lj{vj) 

7 = 1 

(3.6) 

with the linear forms Lj(-) acting on the scalar-valued functions 



L,(v) = J a {/i» + + L Kods 

(3.7) 


In the above formulas, afj, bjj, c u , /,, gf, d i: , hj are sufficiently regular functions defined on 
Cl and on dCl, respectively. 
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Numerous examples fall into the category of problems described by the abstract formu- 
lation (3.1). To mention a few: linear elliptic problems (both single equations and systems), 
linear problems resulting from a one-step approximation in time for evolution problems, 
linear steps of a nonlinear problem solution, etc. 

In this formulation, both the solution u and the test functions v are members of the 
same space X. Non-homogeneous essential boundary conditions are handled by means of a 
standard penalty approach. 


3.2 Finite Element Approximation 

We assume that the domain fl can be represented as a union of finite elements K e , e = 
1,. . . ,M. More precisely 

M 

n= \jK e 

c= 1 

and 

int K t n int Kj = 0 for e ^ / 


Each of the elements K has a corresponding finite dimensional space of shape functions, 
denoted Xh(I\ ), for instance the space of polynomials of order p. The global finite element 
space X h consists of functions which, when restricted to element I< , belong to the local 
space of shape functions Xh(K). Thus the global approximation is constructed by patching 
together the local shape functions in the usual way. 

We shall adopt the fundamental requirement that the global approximation must be con - 
tinuous. As we will see, this requirement leads to the notion of constrained approximations. 
Formally, the continuity assumption guarantees that the finite element space Xh is a sub- 
space of H*(Sl) and, with some additional assumptions if necessary, also a subspace of X. 
The approximate problem is easily obtained from (3.1) by substituting for u and v their 
approximations Uh and Vh- 


Find Uh € Xh such that 

B h {u h , v h ) = Lh{v h ) Vt? A € X h 


(3.8) 


Here 


X h = X h x • • • x Xh{m times) 


(3.9) 


which indicates that the same approximation has been applied to every component of u. 

and L h (-) denote approximations to the original bilinear and linear forms resulting 
from numerical integration. 
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3.3 Adaptivity 

A flowchart of a typical Adaptive Finite Element Method (AFEM) is shown in Fig. 3.1. 
The method consists of first generating an initial mesh and solving for the corresponding 
FEM approximate solution. Next, the error is estimated in some way and based on this 
(usually crude) approximation, one adapts the mesh, i.e., adds new degrees of freedom. The 
approximate problem for the new' mesh is solved again and the whole procedure continues 
until certain error tolerances are met. Obviously, such a procedure requires an estimate of 
the error over each element and a strategy to reduce the error by proper changes in the mesh 
parameters, h and p. 

In our adaptive FEM the new degrees of freedom can be added in two different ways: 
elements may be locally refined or their spaces of shape functions may be enriched by incor- 
porating new shape functions. As noted earlier, in the case of polynomials, this is done by 
increasing locally the degree of polynomials used to construct the shape functions, the first 
case being an /(-refinement, and the second case a p-refinement. A combination of both is an 
(adaptive) h-p FEM. We remark that the process of increasing the local polynomial degrees 
for a fixed mesh size is mathematically akin to increasing the spectral order of the approxi- 
mation and that, therefore, we also refer to h-p methods as “adaptive spectral-element” or 
“A-spectral” methods. 


3.4 Regular and Irregular Meshes 

As the result of local /(-refinements, irregular meshes are introduced. Recall (see [26]) that a 
node is called ixgular if it constitutes a vertex for each of the neighboring elements; otherwise 
it is irregular. If all nodes in a mesh are regular, then the mesh itself is said to be regular. 
In the context of two-dimensional meshes, the maximum number of irregular nodes on an 
element side is referred to as the index of irregularity. Meshes with an index of irregularity 
equal one are called 1-irregular meshes. The notion can be easily generalized to the three- 
dimensional case. (See [7] and literature cited therein for additional references.) 

In the present work, we accept only 1-irregular meshes. In the two-dimensional context, 
this translates into the requirement that a “large” neighbor of an element may have no more 
than two “small” neighbors on a side; in the three-dimensional case, the number of neighbors 
sharing a side may go up to four, while the number of neighbors sharing an edge can be no 
more than two. This is frequently called the “two-to-one” rule (cf. [7]). Examples of regular 
and irregular meshes are shown in Fig. 3.2. There are several practical and theoretical 
reasons to accept only 1-irregular meshes, especially in the context of h-p methods. For a 
detailed argument, we refer to [4]. 

Our restriction to 1 -irregular meshes imposes a simple restriction on the way any h- 
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Figure 3.1: Typical flowchart of an adaptive method. 
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refinement can proceed: before an element is refined, a check for “larger” neighbors must be 
made. If any such neighbors exist, they must be refined first and only then can the element 

in question be refined. 

3.5 Basic Assumptions 

As indicated in the previous sections, the construction of an h-p FEM is based on the 
following assumptions: 

• only 1-irregular meshes are accepted for all /i-refinements; 

• the local order of approximations may differ in each element; 

• the approximation must be continuous. 


3.6 Definition of an Element 

The classical definition of an element is a triple 


{A', A, *>,,*' = 1 — ,A r ) 

where 1\ is the domain of the element (a subset of E 2 or E 3 ), A' is an TV-dimensional space 
of shape functions: 

A B Vi- 1\ ->E, i = l,...,A r 

and Vi, * = 1 A is a set of degrees of freedom, i.e., a set of linearly independent linear 

functionals on X. 

The element base shape functions are understood as a dual basis to Vi- 


y, E X such that 

(Vi, \j) = *> j 1 , • • • , A 


Following this classical construction we define a two-dimensional quadrilateral element as 
follows. In the first step we introduce a two-dimensional master element 

{A, A,<£,,i = 1,..., A} 

The domain K is a unit square, K = [-1, l] 2 . The space of shape functions A is a subset 
of C? P (A), i.e., polynomials of the order p in each variable. We define this subset in such a 
way that A has the following properties: 
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i) Each function in A' can be associated with one of nine nodes of the element: four 
corners, four midside nodes and the centroid; the maximum order of a shape function 
associated with a given node is viewed as the order of this node. 

ii) The base shape functions constituting X will be so-called hierarchical shape functions, 
which means that enriching the element from the order p to p + 1 consists of adding 
some higher order functions to X without modifying the functions already belonging to 

X. 


Space X with these properties is constructed as follows: First we introduce one-dimensional 
hierarchical shape functions on [-1,-1]: 


Xo = |(1-*) 

(3.10) 

xi = ^0+*) 

(3.11) 

X2 = X 2 -1 

(3.12) 

X 3 = X 3 - X 

(3.13) 

\4 = X 4 ~l 

(3.14) 


(3.15) 


(Figure 3.3) The corresponding degrees of freedom can be associated with the two endpoints 
— 1 and 1 and the midpoint 0: 


(<r>0,u) 

= U(-1) 

(3.16) 


= «(1) 

(3.17) 



(3.18) 


where A, are scaling factors. 

Note that the linear functions \o» Xi assume values 0, 1 at the endpoints while all the 
higher order functions vanish at ±1. 

Then for a two-dimensional element we associate the following functions with the sub- 
sequent nodes: 
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Figure 3.3: One-dimensional hierarchical master elements. Degrees of freedom and corre- 
sponding shape functions. 
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i) The bilinear functions 

xo(x)xo(y), xo(x)xi (y), xi(*)xo(y), xiOOxi(y) 

with corner nodes x,y = ±1. 

ii) The functions linear in one direction and higher order in the other, with midside nodes: 

Xa(*)Xo(»)..-.,Xpi(*)Xo(y), 

X2(*)Xi(y)i---iXps(*)Xi(y) 

xo(a-)x2(y),---,xo(x)x P< (y) 

xi(*)x2(y)i---»xi(*)x»(y) 

• The functions at least quadratic in both directions, with the centroid: 

X.'(*)Xj(y)» *1 j = 2 ,...,p 5 

where pi, . . . ,p 4 denote degrees of approximation of midside nodes, p 5 the degree of a 
centroid node (see Fig. 3.4). 

The above set of shape functions is dual to degrees of freedom which are tensor products of 
one-dimensional degrees of freedom given by (3.18). These degrees of freedom can be listed 
as follows: 

• function values at four vertices: 


u( — 1,-1), u(l,—l), u(l,l),u(— 1,1) 


(3.19) 


• tangential derivatives (up to a multiplicative constant) up to p-th order associated with 
the midpoints of the four edges: 



k = 2, . 

• • iP i 

V0(M) 

II 

JO 

• • i Pi 


04 

II 

-se 

••,P3 


k = 2, . 

..,p 4 


(3.20) 
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• mixed order derivatives associated with the central node 

A ‘ A ' &w (0,0) 7 = 2 — p* < 3 - 2i > 

Having constructed the master element we define a subparametric element K. It is obtained 
by mapping of the master element I\ into an actual computational domain: 

K = T (K) 

The mapping T is defined as 

T(x,y) = J2 a 'Ni(x,y) 

i-l 

where A r t , i = 1, . . . , 9 are the regular (Lagrangian) shape functions for the 9-node biquadratic 
element, a t are the desired positions of nodes of K in the computational domain. The space 
of shape functions of K is taken as a space of compositions of T and xp t € -V: 

The degrees of freedom of I\ are defined using s: 

(sPii^j) = f 

with if'} = T 

The definition of a two-dimensional element given above can easily be generalized to the 
three-dimensional case. For completeness, we outline the major steps of this construction: 

We introduce cube subparametric elements: actual elements are images of a cube master 
element K = [— 1, l] 3 under the mapping T : K —* K CR 3 

*‘ = f (3.22) 

j=i 

where N, are second order Lagrangian polynomials associated with 27 element nodes: 8 
corners, 12 midpoints of edges, 6 centers of walls and one central node. 

We equip a master element with the space Xh(I\) of p-th order hierarchical shape func- 
tions defined as a triple tensor product of a set of one-dimensional hierarchical shape func- 
tions \,(') on an interval [—1,1]. Actual shape functions of I\ , A* (A') are as usual compo- 
sitions of the mapping T -1 and shape functions on K: 

X h (K) = {u = uT- 1 | u € X h (K)} . (3.23) 

Degrees of freedom ip associated with the hierarchical base shape functions on a master 
element are defined as follows: 
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• function values at 8 corners: 


u(±l,±l,±l), 


• Tangential derivatives (up to a multiplicative constants) up to p-th order associated 
with the midpoints of edges: 

,d k u 

'ds k ' k = 

where s is a coordinate parallel to the edge, 

• mixed order derivatives associated with centers of walls: 


...... 0 W '“ 

* A ' e7d? 

where s, r is a pair of coordinates parallel to the wall, 
• mixed order derivatives associated with a centroid: 

Qk+l+m 

A -1 A -1 A _1 — — — 

A * A ' Am dx k di ,‘dz”' 

Degrees of freedom of the actual element I\ are defined as 

<<p,u) d = (<p,u) 


3.7 Continuity for Regular Meshes 

One of the fundamental advantages of using the hierarchical shape functions is the ease with 
which they allow one to construct a continuous approximation with locally variable order of 
approximation. A typical situation is illustrated in Fig. 3.5. If elements I\\ and A 2 are 
support polynomials of degree, say, one and three, respectively, then there are at least two 
ways to enforce continuity across the interelement boundary. One way is to add two extra 
shape functions of second and third order corresponding to the middle node A of element K j. 
Alternatively, the same two shape functions may be deleted from element I< 2 - In all these 
cases, a common order of approximation along the interelement boundary can be enforced 
by simply adding or deleting the respective shape functions from the neighboring elements. 
While any of these choices can be made, the results described here employ the “maximum 
rule” in which the higher-order approximation dominates lower orders. Thus, if an element 
is p- refined, i.e., a higher order approximation of degree p is introduced, the neighbors of 
lower order are enriched by the addition of extra shape functions of degree p necessary to 
guarantee continuity of the approximation. 
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Figure 3.5: Continuity by hierarchical shape functions, (cf [9]) 
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3.8 Continuity for 1-Irregular Meshes. Constrained Approxima- 
tion 


Continuity of the approximation on 1-irregular meshes is a more complicated issue. It leads 
to the notion of a constrained approximation that we introduce in this section. 

Consider two adjacent square elements in a 1 -irregular mesh, one having an irregular 
(“hanging”) corner node on the side of the other (Fig. 3.6). The first condition that must 
be satisfied to make the approximation continuous across the common boundary T/ is that 
the spaces of shape functions of elements K e and A'/, if restricted to T/, be identical: 

A'„(A' e )| r , = X h (I<j)\r, (3-24) 

This means, of course, that the orders of approximation for nodes 1 and 2 (Fig. 3.6) must 
be equal. Condition (3.17) is exactly the same as that for regular meshes considered before. 

Denote by *l\,K,( x )i *t>i,K } ( x ) the base shape functions of elements K e and Kj which do 
not identically vanish along T/. 

Assume that the common order of nodes 1 and 2 is p. Since spaces A r / l (A e )|r / and 
X h {I\j)\r } are identical there must exist a unique linear invertible relation between functions 
constituting the basis of these spaces: 

l = 0 ’- -’P (3.25) 

j=0 

In this formula, the functions involved are assigned indices from 0 to p corresponding to 
orders of the functions. Coefficients d R, : are calculated in one of the next sections. The 
superscript d distinguishes between the two generic situations: I\f may be attached to the 
left or the right part of I\ t . Formula (3.25) implies that the degrees of freedom of K t and 1<j 
corresponding to functions involved in (3.25) are also related by a similar equation: take any 
continuous function u>,(a:) such that Uh\h\ € A'/,(A' e ), Uh\K t € Xh(Kj)- Then its restriction 
to Tj can be written as: 


r r 

Wfclr, = Hf ; .V’„K e (*)|r / = 5 Zu, 0 „a r/aOIr, 


(3.26) 


where 

C, = (‘r’i.A'e 1 U h |/\'e) i U i = ^pi,K j i U h\K j') (3.27) 

are values of degrees of freedom obtained for elements I\ e and I< } and v’i.K, are the 

degrees of freedom understood as linear functionals on Xh{K t ) and Xh(I\j). Introducing 
(3.25) into (3.26), we obtain: 

dR H^Kj( x ) lr, = 53 u ^iX-,( x )\rj- (3-28) 

i=Oj=0 i=0 
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Figure 3.6: Two adjacent elements in a 1-irregular mesh. 


49 


Since representation of any function from V p (Tj) in terms of the basis t/’j, K f ( x ) Ir^i j — 
0, . . . ,p is unique, we finally find that: 

u, = '£U/R 1 ,,i = 0 p (3.29) 

i=o 


or, in view of (3.27), 


— Rjii u h\l<c) * 

3=0 


(3.30) 


Formula (3.30) expresses the relation between degrees of freedom of K e and Kj necessary and 
sufficient for continuity of approximation along T/. Necessity was shown above. Sufficiency 
follows from (3.25) and (3.30): 

For any u h such that u A | A ' t € X h (I\ e ), u h \ Kj € X k [K f ), not necessarily continuous, we 
have 


( u /i I ) I r y — y~!( ( r > i.A'« i u h 1a~, ) 0i./v. 1^/ 


is 0 
P 

I 

t=0 

p 

l 
1=0 

V P 


, Uh I Ke ) £ d Rij'Pj,Kj l r / 

3=0 


(UaIA;)^/ = Yl(<PiJx r U h \ K ,)^iJij\Pj - 


= u h \ Kc )^ Kt \Tj 

is 0 J=0 


i.e., (tifc|A«)|r/ = (tlfcUJir/ which means that u/, is continuous along T/. Relation (3.30) 
is referred to as equations of constraints. 

The main idea of a constrained approximation is to replace the degrees of freedom of Kj 
involved in (3.30) by a new set (p iJKr i = 0, . . . ,p related to the old by the matrix d R i: : 

'•PhKj = ]C (3.31) 

«=o 


Then the continuity condition (3.30) in terms of becomes: 

('PiJCjyUklK,) = (*PjJC,iUh\K,) ( 3 - 32 ) 

i.e., the condition which is formally the same as the condition for continuous approximation 
on regular meshes. As a consequence, an element Kj equipped with degrees of freedom 
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can be treated in all procedures of finite element code involving continuity (like assembling 
stiffness matrix) exactly the same way as elements of a regular mesh. 

The degrees of freedom $j,K) defined by (3.31) have a simple interpretation. To compute 
for any ip € Xh(Kj) it suffices to find any continuous extension u/,(:c) such that 
Uh\K, = 4’, u/»|/v, € Xh(I< e ) for which (3.32) holds, i.e.: 

(&.a t r 4>) = (<?>.*«, ( 3 -33) 


In fact, since the action of <Pj,Kt involves only values of along the side r e , we need only 
extend V’ to T e and such an extension, since it must be a polynomial, is unique: 


\p\pj = w(s) = polynomial of s, 


extension of ip to T e *==' w(s) 


(3.34) 


and moreover 


[ ic(s = ±1) , j = 0, 1 


($i,K r 4>) = u h\h t ) = { 


A 7 1 ^(°)’ j = 2 p 


(3.35) 


Illustration of w(s) and degrees of freedom $j,K f is given in Fig. 3.7a. 

We observe that an additional advantage of this approach is that can be associated 

with certain points: endpoints and a midpoint of 1%, i.e., points coinciding with nodes of the 
anticipated neighbor. We call these nodes active (actual) nodes of an element while nodes 
corresponding to original constrained nodes (Fig. 3.7b). 

In case an element is adjacent to two larger neighbors, the transformation of degrees of 
freedom analogous to (3.31) should be performed for the other side of the element as well. 


Interpretation of <Pj.K f shown in Fig. 3.7a suggests a convenient and intuitive way of 
treating a constrained element as if its domain included not only the original square but 
also the adjacent sides of bigger neighbors. This idea gives also a very easy interpretation 
of assembling constrained elements into a mesh: we connect appropriate actual nodes of 
elements as for regular meshes (Fig. 3.8). 

We conclude this section with a discussion of the problem of evaluating the shape func- 
tions dual to the new set of degrees of freedom. Let us write the transformation of degrees of 
freedom resulting from applying formula (3.31) to one or two sides of a constrained element 
in the general form: 

r v>i = H VjRji 1 * € N c 

< (3.36) 

k = 'pi i i 6 IV 
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O 1 — constrained nodes 

• mm* actual nodes 

Node A is both constrained 
and actual 


Figure 3.7: A constrained element: (a) illustration of te(s) and degrees of freedom ( b ) 

active (actual) and constrained nodes. 
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Figure 3.8: Assembling constrained elements into a mesh. 
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where N e is a set of indices denoting degrees of freedom associated with constrained sides 
(i.e., appearing on the left hand side of (3.31) applied to one or two constrained sides of the 
element, whichever case is considered), Af°- in dices of the remaining degrees of freedom, I(i)- 
indices of new degrees of freedom <J5j contributing to a given i Pi, i G A rc . (3.36)2 expresses 
the fact that ^>,-’s not involved in constraints remain unchanged. Relation (3.36) is invertible 
since (3.31) is invertible. The change of basis in any linear space is given by: 

TV 

'fi = Vj^ji 
j= 1 

(3.37) 

Vi = 

j=i 

This implies a change of dual basis according to: 



N 


J = 1 




(3.38) 


Formulas (3.36) correspond to (3.37)!. This means that expressions for the new base shape 
functions can be found by simply transposing the linear relation (3.36) defining 

Transposing (3.36) to obtain results in the following formula: 



^2 RijVj - * 6 A rc , 

i€S(.) 

t/>, , i € N a . 


(3.39) 


where 5(i) = {j G N c :$i contributes to ipj > n (3.36), i.e., t € I(j)}- 

As a simple example of constraints and application of presented formulas consider a 
bilinear element with two constrained sides, i.e., adjacent to two larger neighbors, Fig. 3.9. 
In this case values of degrees of freedom are just values of shape functions at the corner 
nodes and conditions for continuity have an obvious form: a value of a shape function at a 
“hanging” node must be equal to the average of the corner values of the shape function of 
the larger neighbor (Fig. 3.10). Hence, recalling the interpretation of <£,•’ s (Fig. 3.7a), we 


54 



write equations (3.36) for element considered as follows: 


<r>l = 

$1 




1 _ 

1 _ 


<r>2 = 

2* 1 + 

2^ 


V>4 = 

1 . 
2^ 


1 _ 

+ 2* 4 

= 



$ 3 


In this case, N c = {1,2,4}, N a = {3}. The matrix form of (3.40) is: 



According to (3.3S)i the base shape functions are given by 



(3.40) 


(3-41) 


(3.42) 


(3.43) 
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Figure 3.11 presents the constrained base shape functions 0,-. Relation (3.43) could also be 
formally obtained using the genera] expression (3.39). 


3.9 Calculation of the Element Load Vector and Stiffness Matrix 


For simplicity, we restrict ourselves to the case of a single equation. In the case of systems 
the same procedure is applied for every linear form (3.7) and bilinear form (3.4). 

The element load vector and stiffness matrix are defined as: 


< 


b, 

Bn 


Lh,K (^i) , 
Bh,k (V\,Vb) i 


(3.44) 


t,j = 1 , . . . , jV, N being the number of degrees of freedom of an element, t/>, base shape 
functions corresponding to actual degrees of freedom 

Since constrained shape functions xpi are related to the usual shape functions xj.\ by a 
linear transformation (3.39), 6, and B, } can be expressed as follows: 

' b, = RijLhAtt'j) for i G N c , 

i€S(t) 

< 

k bi = LkjcWi) for i G N a , 

r Bij = E E RikRitBkMk^i) if 

kesu)tesu) 

(3.45) 

B,j = E R>kBh,i\{ipk,il>j) if 1 € N c ,j 6 N a , 

t€S(.) 

i 

= E R ikBh, kWM if t G N a ,j € N c , 

S(j) 

„ Bij = Bk,s(tl>i,rl'j) if i,j € N a 

where Lh,k, Bh,i< are restrictions of forms (3.7) and (3.4) to the element I \ , sets of indices 
N c , N a , S(i) were defined in (3.36) and (3.39). 
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Figure 3.10: Conditions for continuity for a constrained bilinear element. 
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3.10 Constraints in the One-Dimensional Case 

In the case of irregular meshes, continuity has to be enforced by means of the constrained 
approximation. To fix ideas, consider the generic, one-dimensional case shown in Fig. 3.12. 
The approximation on the small elements [—1,0] and [0,1] must match the approximation 
on the large element [—1,1]. 

We first choose the scaling factors A p in (3.18) in such a way that the corresponding 
shape functions for the one-dimensional master element have the following form 

Xo = i(l-0 

Xi = ^(1+0 ► (3.46) 

e-1 P = 2,4,6,... 

X P (0 = < 

£ p — £ p = 3,5, 7, . . . , 

Assume next that all degrees of freedom for the large element are active. The question 
is: what degrees of freedom must the small elements take on in order that the functions 
supported on the two small elements exactly coincide with shape functions of the large 
element? 

From the fact that (3.18) is a dual basis to (3.46), we get 

V>,(X>) = T P'- *= 1 P = 2,3, . . . (3.47) 

Ap 

and therefore A p = p!. 

The transformation map from [—1, 1] onto [—1,0] is of the form 

x = -\ + \t ( 3 ' 4S ) 

with inverse £ = 2x + 1. This yields the following formulas for the shape functions 
e Xp<P — 0,1,2,... For the (left-hand side) element [-1,0] (recall definition of a subpara- 
metric element). 

*Xo(*) = ~ x 

( X i{x) = x + 1 (3.49) 

e Xp (x) = 1 — (2x + l) p p = 2,4,6,... 

' Xp (:r) = (2:r + l) p — (2*+l) P = 3,5,7,... 

and the corresponding formulas for the degrees of freedom are 

<'<p 0 ,u> = u(-l) 

<'<pj,u> = u(0) 
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Figure 3.12: Derivation of the constraint coefficients. 



, 1 d?u ( \\ 

< >~ 2 P P \dxt> \ 2 ) P ’ 


(3.50) 


Now let u(z) (x = £ for the master element) be any function defined by the shape 
functions on [—1, 1], i.e., 

u(x) = '52Vi( u )Xi( x ) ( 3-51 ) 

9=0 

In order to represent u(x) for x £ [—1,0] in terms of the shape functions on [—1,0], we 
have to calculate the values of the degrees of freedom (3.50). We get 

k 

< e <?o ,u> = spo(ti) <* <*,Xo >+!>,(«) <' 


where 


For p > 2, 


where 


= < <Po, u > 

< l ,v > = <,?<>(«) < l ‘PlfXo > +<r’l( u ) ‘r’lfXi > 

k 

+ ^2^p q (u) <‘ <puX, > 

9=2 

1 

= - < v>,,u > + £ l R ql < > 

1 7=2 


Rql — ^ ^ 

1 if 9 is even 

0 otherwise 


< / U > = ^o(^) <* XO > + H ¥>pi X? > 

9=1 

k 

= 0 + 51 ‘Rqp <<fq,U> 

9=1 


Rqp — ^ ^PpiXq > 


for 9 < p 


= I,_ ir nM = — for ,>p 

1 2< ( ' W 2 * f>H« - J>)! 

The same procedure applied to the right-hand side element [0,1] yields the following: 
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The tratisformalion from [—1,1] onto [0,1] 


with inverse £ = 2x — 1. 

The shape functions r \ p ,p = 0,1,2,...: 


(3.55) 


r \'o(*) = 1 “ * 

(3.56) 

r X,(x) = x 

(3.57) 

_ j 1 - (2x- l) p peven 

* p 1 | (2x - l) p - (2x - 1) p odd 

(3.58) 


The degrees of freedom 'V P : 


< Vo, u> = 

u(0) 

< Vi , u > = 

u(l) 

r 1 

d p u /1\ 

^•" >= 2 y 

dx^ \2/ 


(3.59) 


The constraints 


where 


1 1 k 

< Vo, a >= - < 9?o,u > +- < 9 j,u > + '52 T R q O < pq ,1 

2 9=2 


{ 1 if q is even 
0 otherwise 


< Vi, u >=< > 

and for p > 2 

k 

< 'r l V' u ' > = ^- ^ Rqp < "~ ^Pqi u ^ 
9=2 


(3.60) 


(3.61) 

(3.62) 

(3.63) 
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where 


— ( Vpi Xq 

[ 0 for q < p 


(3.64) 


' R 


QP 


I/V 

Oq 


for q > p 


(3.65) 


Finally, the conditions for matching the approximation on one-dimensional master ele- 
ment [—1,1] and on the two small elements [—1,0], [0, 1] can be rewritten as: 


Vi = i = 0,...,p 


(3.66) 


j=0 


where the superscript d — l (left) or r (right) distinguishes between the two small elements. 
Formula (3.66) relating degrees of freedom of the two elements is equivalent to a correspond- 
ing relation between the base shape functions d \i and dual to V, and 9,: 


Xi = tl dR H d X J, 

j - o 

Arrays '/? 9P and r /? 9P , q,p = 0, ... ,5 are presented in Fig. 3.13. 


(3.67) 


3.11 Constraints for Two-Dimensional Subparametric Elements 

Since the shape functions for the 2-D master element are defined as tensor products of 
the 1-D functions, the results for the 1-D case hold exactly in the same form in the 2- 
D situation, the only difference being that the calculated constraint equations have to be 
applied to the proper degrees of freedom (see 3.31). It follows from the definition of the 
subparametric elements that the constraints coefficients are exactly the same , even w'hen 
the elements have curved boundaries. This follows from the fact that the shape functions 
behavior in a subparametric element on a part of its boundary depends exclusively upon the 
deformation of the part of the boundary, and therefore, any relation defined for the shape 
functions in the generic situation carries over immediately to the case of two small elements 
sharing an edge with a large element, so long as the deformation of the edge is identical in 
all three elements. The situation is illustrated in Fig. 3.14. 


3.12 Constrained Approximation in a Three-Dimensional Case 

A three-dimensional constrained approximation exploits basically the same concepts as we 
developed for the two-dimensional case. The constraints, however, in this case have a more 
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Figure 3.13: The constraint coefficients for a sixth order approximation. The unfilled coeffi- 
cients are zero. 
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complex form. Moreover, there are two different types of constraints: one resulting from 
constraining only an edge of the element and the second resulting from constraining the 
whole wall of the element. 

(a) an element may have a larger neighbor across the wall, Fig. 3.15a 

(b) the neighbors across walls are of the same size as the element considered, but the 
neighbor across an edge is larger, Fig. 3.15b 

These two situations must be considered separately. 


3.13 Constraints for a Wall 


We construct equations of constraints and transformation to new degrees of freedom and 
base shape functions in essentially the same way as we did in the two-dimensional case. 

We consider two adjacent cube elements I( e and Kj of a 1-irregular mesh, one of them, 
A'/, two times smaller than I< t . The wall Wj of Kj is attached to one quarter of the wall 
W e of the larger neighbor, Fig. 3.16. The first condition for continuity of approximation is 
that spaces of shape functions Xh(I\ e ) and Xh(Kj) if restricted to Wj must be identical: 

X h {K } )\w, = A' fc (A' e )| HV (3.68) 


and it is satisfied since we accepted the “maximum rule”. 

In the reasoning below we use two-indices notation for shape functions of Kj and K e 
which do not vanish identically on Wj : an d V’O.K'/ - Restrictions of these functions to 

Wj are of the following form (Fig. 3.15) 


i wj = Xt (* ) • \] (y ) *, j = o, • • • * p 

vSj./viiv, = di Xi{0- dl Xi(v) = 


(3.69) 


where x, and dl \,, dl \, are one-dimensional hierarchical shape functions considered in the 
section on “Constraints in One-Dimensional Case,” d\, d 2 indicate which of the two generic 
situations considered there should be applied. 

Functions \ Pij,K,\W] constitute two bases of spaces Xh{Kj)\w, = Xh{K e )\\Vj, 

therefore there must exist a unique linear invertible relation between them. Using the trans- 
formation (3.67) between x. W an d dl x(0i an ^ the same for the y-direction we easily find 
that this relation is of the form: 

*<*.!••, -££**» ** jl j\W j (3.70) 

it=0 (=0 
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b) 



Figure 3.15: Three-dimensional constrained elements: (a) an element with constrained walls; 
(b) an element with a constrained edge. 
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K, 


Figure 3.16: Two adjacent three-dimensional elements ina 1-irregular mesh having a common 
wall. 
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which is analogous to formula (3.25) obtained for two-dimensional case. 

By exactly the same reasoning as in that case we can show that introducing new degrees 
of freedom p,j,K } by 

*PH UK, = ( 3 ‘ 71 ) 

k = 0 (=0 

we obtain the necessary and sufficient condition for continuity of approximation along the 
wall of the form: 

{$ij,K } ,Uh\K } ) = i\Ke) ( 3 - 72 ) 

for every Uh such that Uh\h\ € A^(A' e ), Uh\!<j £ Xh(Ej)- 


3.14 Constraints for an Edge 


Consider now a situation when a smaller element I<j has a larger neighbor across the edge 
Ej even though the neighbors across the walls adjacent to Ej are of the same size as A/, 
Fig. 3.17. In this case only the shape functions associated with the edge Ej are involved in 
constraints. 

Let i/’i h\ i 4’i,Kj be the base shape functions of I\ e and Kj which do not identically vanish 
on Ej. Again, continuity of approximation requires that: 

A'»(A'.)Ie, =A'»(A'/)|e,. (3.73) 

The two sets of functions V’l.K.lEp are the basis of the above spaces, and since. 

{ ^i,K t \E, = X.(z) 

(3.74) 

= d Xi(0 

and \,-(x), d \i(0 are related by (3.67), we find that 

i\K'\E, = '£ d Rii r l’j,Kj\E } - (3.75) 

J=0 


By the arguments used before we introduce a new set of degrees of freedom pi.Kj such 
that: p 

pi.K, = Y dR > i ( 3,7S ) 


j=0 


and as a consequence we obtain the necessary and sufficient condition for continuity along 

Ej: 

(p « ,K j i ^ ft|/\ /) {'r’i.A'* i 

for all Uh such that Uh\K c € A'/i(A' e ), Uh\K f € A^(Ay). 


(3.77) 
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Figure 3.17: Two three-dimensional elements in a 1-irregular mesh with a common edge. 
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3.15 The Constrained Base Shape Functions 

We collect transformations (3.71) and (3.70) obtained for all constrained walls of the element 
and all constrained edges not belonging to constrained walls. As a result we end up with 
the relation which, after introducing a uniform one-index numbering of degrees of freedom, 
can be written in a general form: 


’ w = Y 'SjRj* 

, ieN c 

;=/(•) 

1 


II 

, ieN a 


where N c is a set of indices denoting degrees of freedom associated with constrained walls 
and edges of the element, ^“-indices of the remaining degrees of freedom, /(*> indices of &’s 
contributing to a given ^ in (3.78). Coefficients are expressed by dl R, : , d3 R, : in a way 
indicated by (3.76) or are just equal to d R tJ if constraints result from (3.78). 

By the arguments used before the base shape functions 4 \ corresponding to are given 
by a linear transformation transpose to (3.78): 


’ 4\ 

= Y, r ‘^'j 

j€S(.) 

, t € A c 

(3.79) 

. * 

= 

, ieN a 



where 5(t) = {j € N c | i € /(.))}• 

3.16 Interpretation of Calculation of the Load Vector and 
Stiffness Matrix 

As in the two-dimensional case, degrees of freedom have a simple interpretation: to 
compute € A^(A') we can uniquely extend the function as a polynomial to 

constrained walls and edges extended such that they coincide with walls and edges of the 
anticipated larger neighbors. Then, appropriate degree of freedom v?i of the larger neighbor 
should be applied for the extension of 4 '. All generic situations are presented in Fig. 3.18. 

The algorithms for transforming the usual element load vector and stiffness matrix to 
those corresponding to s are exactly the same as in the two-dimensional case. The only 
difference is that they use more complex formulas (3.79). 

Also, the arguments that we used to extend the concept of constrained approximation to 
two-dimensional subparametric elements apply in the three-dimensional case. 
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3.17 Concluding Remarks on Constrained Approximation 

• A formal definition of an element states that it is a triple: 

{A', A r /i(/\ ), {<^i},=i n ) i ( 3 - 8 °) 

where K C J? 2 (J? 3 ) is a domain of the element, Xh{X) the space of shape functions, 
{s?i}»=i N a set of linear functionals on Xh{ K) called degrees of freedom. 

From this point of view the approach that we propose for enforcing continuity on 1- 
irregular meshes is equivalent to constructing a new element: we define new degrees of 
freedom. 

• A finite element code using constrained elements must include only three non- 
traditional algorithms: 

— a procedure identifying kinds of constraints for a given element, 

— an algorithm transforming the usual load vector and stiffness matrix to those 
corresponding to actual degrees of freedom £>,, 

— the procedure transforming the finite element solution in terms of ipi s (obtained 
from the solver) to values of usual degrees of freedom, i.e., the procedure perform- 
ing the calculations indicated by (3.36) or (3.78). 

These three algorithms involve complex logical operations; however, once they are coded, 
they may be used as “black boxes” by a user not familiar with their content. The rest of the 
code is unaffected by the constrained approximation and therefore it may be developed in a 
standard way. 

3.18 Some Details Concerning the Data Structure 

In the classical finite element method, elements as well as nodes are usually numbered consec- 
utively in an attempt to produce a minimal band within the global stiffness matrix. When 
the program identifies an element to process its contribution to the global matrices, the 
minimal information needed is the node numbers associated with the element. Adaptive re- 
finement and unrefinement algorithms require much more information on the mesh structure 
than the classical assembly process. 

First of all, we introduce the notion of a family. Whenever an element is refined a new 
family is created. The original element is called the father of the family and the four new 
elements are called its sons. Graphically, the geneology on families can be presented in a 
family tree structure as illustrated in Fig. 3.19. 
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An examination of refinement and unrefinement algorithms (see [7] for details) reveals 
that for a given element NEL, one must have access to the following information: 

• element node numbers 

• element neighbors 

• the three structure information, including: 

— number of the element family 

— number of the father 

- numbers of the sons 

- refinement level (number of generation) 

For a given NODE we also require, 

• node coordinates 

• values of the degrees of freedom associated with the node 

In general, some information is stored explicitly in a data base consisting of a number of 
arrays, some other information is recovered from the data base by means of simple algorithms. 
A careful balance should be maintained between the amount of information stored (storage 
requirements) and recovered (time). 

The following is a short list of arrays used in the data base: 

1. The tree structure is stored in a condensed, family-like fashion [26], [7] in two arrays 

NSON(NRELEI) 

NTREE(5,MAXNRFAM) 

where NRELEI is the number of elements in the initial mesh and MAXNRFAM is the 
anticipated maximum number of families. For an element NEL of the initial mesh, 
NSON(NEL) contains its first son number (if there is any). For a family NFAM, 
NTREE(1,NFAM) contains the number of the father of the family while the other four 
entries NTREE(2:5, NFAM) are reserved for the “ first-born ” sons of the sons of the 
family (the first-born “grandsons” of the father). 

2. The initial mesh neighbor information is stored explicitly in array 

NEIG(4, NRELEI) 
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Figure 3.19: A tree structure and the natural order of elements: 4, 5, 12, 13, 14, 16, 17, 18, 
19, 7, 8, 9, 10, 20, 21, 22, 23, 3. (after [9]) 
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containing up to four neighbors for each element of the initial mesh (elements adjacent 
to the boundary may have less neighbors). 

3. For every active element, up to nine nicknames are stored in array 

N0DES(9,MAXNRELEM) 

where MAXNRELEM is the anticipated maximum number of elements. 

For a regular node, the nickname is defined as 

NODE* 100 + NORDER 

where NODE is the node number and NORDER the order of approximation associated 
with the node. 

For an irregular node, the nickname is defined as 
NORDER 

where NORDER is again the order of approximation corresponding to the node. 

4. For a particular component IEL of a vector-valued solution, the corresponding degrees 
of freedom are stored sequentially in array 

U(MAXNRDOF,IEL) 

where MAXNRDOF is the anticipated maximal number of degrees of freedom. Two 
extra integer arrays are introduced to handle the information stored in array U. Array 

NADRES(MAXNRNODE) 

contains for every node, NODE, the address of the first from the degrees of freedom 
corresponding to NODE in array U. If K = NADRES(NODE) is such an address, the 
address for the next degree of freedom can be found in 

NU(I\) 

and so on, until NU(K)=0, which means that the last degree, of freedom for a node has 
been found. The parameter MAXNRNODE above is the anticipated maximal number 
of nodes. 

5. The node coordinates are stored in array XNODE 

XNODE(2, MAXNRNODE) 
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The rest of the necessary information is reconstructed from the data structure by means of 
simple algorithms. These include: 

— calculation of up to eight neighbors for an element 

— calculation of local coordinates of nine nodes for an element determining its ge- 
ometry (the irregular nodes coordinates have to be reconstructed by interpolating 
regular nodes coordinates) 

— recovery of the tree-structure related information, e.g., level of refinement, the 
sons’ numbers, etc. 

— an algorithm establishing the natural order of elements 

During the h and p refinements and unrefinements, both elements and nodes are created 
and deleted in a rather random way. This makes it impossible to denumerate them in a 
consecutive way, according to their numbers (for instance, as a result of unrefinements some 
numbers may be simply missing). Thus a new ordering of elements has to be introduced 
which is based on some scheme other than an element numbers criterion. In the algorithms 
discussed here, we use “the natural order of elements' " based on the initial mesh elements 
ordering and the tree structure. The concept is illustrated in Fig. 3.19. One has to basically 
follow the tree of elements obeying the order of elements in the initial mesh and the order 
of sons in a family. 

The natural order of elements may serve as a basis for defining an order for nodes and, 
consequently, for degrees of freedom, when necessary. 

For a detailed discussion of the data structure as well as a critical review of different data 
structures in context of different /i-refinement techniques, we refer again to [7]. 
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4 Adaptivity 


The main advantage of an h-p finite element method is the possibility of adapting the mesh 
to features of the approximated solution. Adaptivity should lead to enriching the current 
mesh only where it is needed, i.e., where the accuracy is not sufficient and where the new 
degrees of freedom cause the best improvement of the solution. 

There are two basic steps in the process of adapting the mesh. The first is the a posteriori 
estimation of errors of the current approximation. This estimation provides the necessary 
local information about the quality of the solution required to drive the adaptive strategy. 
In the second step, based on the knowledge of the errors, we refine the mesh: break or enrich 
the elements. The rules for making the decision as to which elements should be refined or 
enriched play the key role in generating optimal meshes. They are usually referred to as 
adaptive strategies. 

In the following, we first present error estimation techniques which have been imple- 
mented for the compressible Navier-Stokes equation. Then we discuss extension of these 
techniques to calculate the directional adaptation indicator. The h-p mesh adaptation 
strategies are described in the last subsection based on these indications. Note, that the 
error estimation and h-p adaptive techniques were developed in the previous year, while the 
directional adaptation indicator is a recent development in this project. 


4.1 Error Estimation Techniques 

In general, there are two major classes of error estimation techniques: interpolation error 
estimates and residual error estimates. The former group of methods, interpolation methods, 
provide a rather inexpensive approach for estimating the numerical error. This approach, 
however, is usually not very accurate and only provides a relative indication of where large 
errors exist. The latter class of methods, residual methods, are typically much more accu- 
rate but are also much more expensive to use. During the course of this project we have 
experimented with both classes of methods. 

4.1.1 Interpolation Error Estimate 

This method of error estimation employs well known a priori estimates of the interpolation 
error of finite element approximations. Such estimates are given by the following formula: 

II" “ u /llo,K = c Nil,* ( 4<1 ) 
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where uj is an interpolant of u. Of course we are not interested in the accuracy of the 
possible interpolation of the exact solution u, but rather in the accuracy of the finite element 
approximation of the problem u h - Still, numerical experiments indicate that the errors given 
by (4.1) can be considered a rough indication of the accuracy of uh and can serve as a basis 
for mesh adaptation. 

The major advantage of this method is that it is computationally inexpensive, problem 
independent, and easy to implement. Yet its rather poor quality is a reason that we use it 
only if other techniques are not available or if we need only a very rough estimate of the 
error. 

4.1.2 Residual Error Estimate 

The idea behind residual a posteriori error estimates can be outlined as follows. We substitute 
the existing finite element approximation u h into the original statement of the problem being 
solved. Since u h is not an exact solution we obtain a certain residual r h which could be 
measured in a suitable global norm. For instance, it could be the norm of the space dual to 
the space containing the solution A : 


INI = sup 

v£A' 


(n,,n) 

IMI 


(4.2) 


for which we are guaranteed that it exists. The exact solution of our original problem is, 
however, usually unavailable and we can only try to estimate the value of this expression. 
The techniques leading to such an approximate evaluation of ||r A || are referred to as residual 
error estimates. They express || || as a sum of element contributions which we call local 
error indicators and they also reflect the local accuracy of the solution (i.e., for each element). 

The element residual method was originally developed for symmetric elliptic bound- 
ary value problems. The method was recently extended to a class of nonsymmetric but 
symmetrizable problems which includes compressible flow problems [24]. In the following 
discussion we will provide details for the nonsymmetric version of the method and give onljr 
a general outline of the method for the symmetric case. For details about estimating ex- 
pression (4.2) by local element contributions, we refer to [25]. The presentation below is 
extracted from [24]. 


Element Residual Method 

Given a domain ft C R N (we assume TV = 2 for notational simplicity) we consider a general 
variational boundary value problem in the form 
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{ 


Find u 6 X such that 
B{u,v) = L(v) for every v € X 


(4.3) 


where 


x = H\n) = H\n) x ... x H'{n) 

' v ' 

n times 


B(u,v) = ^2 Bijiu^vj) 

*J = 1 

l(v) = £i>,) 

j=l 


(4.4) 


with the bilinear forms B tJ and linear forms Tj defined as (omitting superscripts for nota- 
tional convenience) 


B(u,v ) = 


+ 


/ / ^ du dv du , dv 1 

akl dx ( dx k ^ k dx k V ^ tU dx e CU J 

r ( du dv \ 

/ <b 3 —v + d s u— + c t uv>ds 

Jan [ ds ds j 


dx 


(4.5) 


L(v) = IJfv + 


f dv\ 


dx 


(4.6) 


+ 


/ /a 

dan 


r ds 


For each pair of indices ij = 1 , . . . , n, a kh b k , d t , c, /, 5 / are functions specified in f) and b„ 
d a , c 4 , / a are functions specified on the boundary dCl. The normal and tangential derivatives 
on the boundary are defined as 


du 

du 


du 

dn 

= ^ 

+ 


du 

du 


du 

ds 


+ 

a - 
dx 2 


(4.7) 


where (n!,n 2 ) are components of the outward normal unit vector n. 
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Systems of type (4.3) include not only classical elliptic equations of second-order but 
also arise naturally as “one time step problems” from different time discretization schemes 
applied to parabolic or hyperbolic equations. The boundary integrals in (4.5) permit the 
implementation of different boundary conditions (including Dirichlet boundary conditions 
via the penalty method). 

Replacing X in (4.3) with a finite dimensional subspace Xh, p of X we arrive at the 
approximate problem 


Find Uh, p € Xh, P such that 


{ B(u htP ,v) = L(v) VveX h , p 


(4.8) 


Indices h and p refer here to the use of an arbitrary h-p adaptive finite element (FE) meshes, 
with locally varying mesh size h and spectral order of approximation p. 

It is our goal to propose and investigate here a general method for estimating the relative 
residual ert'or corresponding to (4.8). More precisely, considering the enriched space X h , P +\ 
corresponding to the same mesh but with local order of approximation uniformly increased 
by one, we define the relative residual error as 


sup 

vzX h lP + 1 


\B{u h , p ,v) - I(t>)| 

IMI 


(4.9) 


The choice of norm ||v|| is unfortunately not unique. Two important special cases are, 
however, of interest: the symmetric case , when B is symmetric and positive definite, and the 
symmetrizable case when B can be made symmetric by an appropriate change of variables. 


Symmetric Case 

When the bilinear form B is symmetric and positive definite and the energy norm 


He = B(v,v) (4-10) 

is selected in (4.9) the residual error is equal to the relative error between u Kp and u/i.h- i, 
the FE solution corresponding to the enriched space and measured in the energy norm. 


sup 

veX h.p+i 


\B(uh , P ,tQ - Lj t>)| 
IMIe 


II Uh,p - ttfc.p+l||E 


(4.11) 
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The principal idea behind the proposed error estimate is to interpret (4.11) as a variational 
formulation of an elliptic problem, transform the bilinear form B into the typical form for 
elliptic problems, and finally apply the element residual method presented in [25]. 

Formally, we proceed as follows: 


Step 1: Transform formulas (4.5) and (4.6) into the typical form for elliptic equations. 

2 du dv 


B(u,v) = 


+ 


+ 


L(v) = 


/ I A au ov 

Jn \^ aU dx ( dx k 

L{ b 'Ts v+d -*Ts + { c - + t/‘ n ) m '} ds 


(4.12) 


f-iM 

h\ ! t,dx,) 


vdx + 


/ ('• + £ 

ha \ 


gtfi{ vds 


Step 2: Apply the element residual method to the modified bilinear and linear forms result- 
ing in the estimate 


\u h , p - u h<p+1 \\ E 


— ^51 IIv’aIIe.a- j 


(4.13) 


where the error indicator function is the solution to the local problem 
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f Find <p K € X° hp+l (I<) such that 


B K {<p K ,v) = 



[ for every v 6 -X’° p+1 (A') 

Here X° hp+1 {I<) is the kernel of the h-p interpolation operator defined on the ele- 
ment enriched space X h tP +\(K) or the so-called space of element bubble functions 
and the element bilinear form B K is defined as the element contribution to (4.10). 
Finally, the symbol [ ] denotes the average flux defined along the interelement 
boundary and evaluated using both the element and the neighboring elements val- 
ues of derivatives and coefficients ajf, (if they are discontinuous). The element energy 
in (4.13) is defined using the element bilinear form Bk- 

Step 3: Integrating by parts transforms the element bilinear form and the right-hand side 
of the local problem into a form consistent with the initial formulas for B and L. 
Thus, we arrive at the following formulas 
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(4.15) 


«,i=i 

Lk{v) = 

i=i 


where 



+ / f,vds 

JdKndQ 

The final form of the local problem is derived as follows: 


f Find <p K € X° hp+1 (I\) such that 


Bk{<p k ,v) = L h {v) - Bfc(uh, p ,v) 


+ 


£ / 

, J8K\ea 


y' a *i^i 


Vjds 


(4.17) 


Nonsymmetric and Symmetrizable Problems 

Formally, formula (4.13) can be used for nonsymmetric problems as well, as long as the local 
element bilinear forms Bk are positive semidefinite, i.e., 

Bk(v>k^k)> 0 ( 4 - 18 ) 
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This happens if the symmetric contributions to Bk dominate the unsymmetric ones (result- 
ing usually from the first-order terms). The global bilinear form B is then automatically 
semipositive and, with the correct boundary conditions, it is positive definite. This guaran- 
tees the well-posedness of the problem. 

Another interesting case is when the bilinear form is nonsymmetric but it is symmetriz- 
able in the sense that a matrix- valued function A.o(x) exists (the so-called symmetrizer 
introduced earlier) such that a new bilinear form B defined as 

B(u,v) = B(u,A 0 v) (4-19) 


is symmetric. 

If, in addition, the symmetrized bilinear form B is positive definite, then the error esti- 
mation technique can be extended to this case as well. 

Introduction of the symmetrizer does not effect the construction and solution of the local 
problems. It only helps identify the norm for the space Xh , P + 1 in (4.9) and affects the 
evaluation of the error estimate. Using the same definition of element bilinear and linear 
forms Bk,Lki we proceed as follows: 

Step 1: Use the orthogonality of the residual to the Xh, P space, 


B{u h>v ,v) - L(v) = B(u hlP ,rp) - 

(4.20) 

where 

ip = V - U hiP v 

where II^p denotes the h-p interpolation operator (see [25]). 

(4.21) 


Step 2: Decompose the bilinear and linear forms according to formulas (4.12) introducing 
the average flux interelement boundary terms 

B{u hiP ,v) - L{ v) 

Step 3: Introduce the solutions to the local problems 

B(uh tP ,v) - L(v) 

k 

where xp K is the restriction of %p to element K. 


£ a ^. nk 
k,i=\ uxi 


ip- 'ds 


(4.22) 


(4.23) 



Step 4: Introduce the symmetrizer and use the Cauchy-Schwartz inequality for the sym- 
metrized form to estimate the error 

B{u hyP ,v) - L(v) - Y^B K (<p K ,A 0 Ao'ip K ) 

1< 


= < Y.^^K^K) lB K{A 0 i ^ K ,AQ^ K )^ (4.24) 


K 

< c 


^Bk{v>kiAqV>k) 

L K 


B{AZ'v,v)2 


Here C = max a - Ck where for every element K , Ck is identified as the norm of 
(/ - IU,) operator with respect to the element energy norm defined as 

IMIe.k = -EMV”,>’) < 4 - 25 ) 


(see [25] for a detailed discussion of C). For undistorted meshes C is close to 
one (independent of the order of approximation) and in practical calculations is 
neglected. 

Identifying the global energy norm for v in (4.9) as the sum of (4.26) we arrive at 
the final estimate of the form 


sup 

V€-X\, p +i 


|jg(u/,, p .i>) - L(i;)| 

IMI 


< 




1 

2 


(4.26) 


Example: Taylov-Galcrkin Method for Euler Equations 


Recall that the variational formulation of the Taylor-Galerkin method is of the general form 
(4.3) with the bilinear form defined as follows. 


B(U U+ \V) 



+ boundary terms 


(4.27) 


The form of boundary terms present in the formula for the bilinear form depends on boundary 
conditions. 

The formulation is nonsymmetric. However, it is known (see [16,14]) that there exists a 
symmetrizer A 0 = A 0 (u), see Fig. 2.1 (Hessian of the entropy function for Euler equations), 
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such that 

1. a 0 = aZ>o 

2. {A 0 A,) T = A 0 A t > 0 


(4.28) 


Based on (4.28), one can easily verify that (with a proper treatment of boundary conditions) 
the bilinear form 

B{U,V) = B(U,A 0 V) (4.29) 

is symmetric, provided the derivatives of the symmetrizer Aq are negligible, i.e., 



(4.30) 


Example: The Momentum Step of the Navier-Sf okes Equations 


The momentum step in the two-step procedure outlined in Section 2 involves solving the 
system of equations: 


mj n+1 - 0&t £ r"* 1 = rv] + (1 - /3)At £ T? jti 
1=1 


(4.31) 


Equations (4.31), if rewritten in terms of the velocity components, reduces to a system of 
two symmetric, elliptic equations. Unfortunately, in order to comply with the conservative 
form of the equations, (4.31) must be solved in momentum components. 

The variational formulation of (4.31) does not result in a symmetric problem but the 
bilinear form may be symmetrized using the symmetrizer 


Aq — 


■ 1 

p 

0 


0 

1 

P . 


as this transforms the problem to a symmetric formulation in velocity components. 


Example: The Energy Step for Navier-Stokes Equations 
The energy step involves solving the equations: 

e " +1 - (trf'ur 1 + «T?' 

•=i \j=i j 

= e" + (1 - /*)A ft |S>”u" + 

t=i V=1 


(4.32) 
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The variational formulation of (4.32) is not symmetric. However, since rewriting (4.32) in 
terms of temperature T results in a symmetric diffusion equation, and since e = c v pT + m 2 / 
2 p, the factor Aq = 1/p is a suitable symmetrizer of the problem. 

The extension of the element residual error estimation to the implicit slash explicit 
method (which is based on one-step Taylor-Galerkin Formulation) is straightforward: first 
we calculate the error indicator function by solving the local problem (4.23) for each element, 
then use (4.26) to compute the error indicator for that element. The bilinear form is obtained 
from the variational formulation of the problem as before, and is of the general form (3). The 
same symmetrizer used for Euler equation can also be applied for Navier-Stokes equations 
(cf Hughes’ paper). It should be noted that, although the algorithm extends in a neutral 
way, the theoretical work for the Navier-Stokes equation is still not complete. 


Numerical Examples 


In this section, two example problems illustrating these techniques of error estimation are 
presented. Note, that these are rather simple examples designed to illustrate the basic 
ideas presented here on relatively worse meshes. More practical applications we presented 
in Section 7. The results take the form of plots of the error estimates and effectivity indices 
as well as global effectivity indices and standard deviations. These quantities are defined as 
follows: 

1K = it (4-33) 

where jk is the effectivity index for element K, 6k is the estimated error and |||e|||A- is 
the actual element error in the coarse mesh approximation (comparing the coarse mesh 
approximation with either the analytic solution or the approximate solution on a mesh 
of uniformly increased polynomial order). Additionally, we introduce a discrete measure 
(weight) u>k defined according to 


u>k = 


\l 


IIHIP 


(4.34) 


With this definition, the global effectivity index becomes 

(?'•) 


7 2 = 


< P 


= 7k u k 


(4.35) 


K 


Now classical statistics suggest a standard deviation a (with respect to the measure) as a 
method to quantify the ability of the estimates to predict an appropriate distribution of 
error. The standard deviation is defined as: 

L 2 


^ = E (•& - -T 2 )** 

K 


(4.36) 
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In order to eliminate any global constants that may be missing from our estimates, we 
normalize the element effectivity indices by dividing them by the global effectivity index: 


'"'■“IMIIk ' 7 

which results in a standard deviation defined according to: 

a 2 = 53 (lK “ O*"* 

K 


(4.37) 


(4.38) 


Example 1: Inviscid Flow Over a Blunt Body 

We used the Taylor-Galerkin method described in Section 2 to model the flow over a blunt 
body with Mach number M = 6. Figure 4.1 shows the density contours of a steady-state 
solution obtained on a uniform mesh of 16 x 16 linear elements. Figures 4.2a and 4.2b present 
distributions of the error indicators 6 K (obtained using (4.26)) and the normalized effectivity 
indices 7 K (4.37). Since the exact solution to the problem is not available, the exact errors 
are not known. For this reason we computed the effectivity indices 7K = ^k/IIMIIk, us > n 6 
instead of the true errors |||e||| K , the errors understood as a difference between the actual 
finite element solution and the solution obtained by performing one time step on the mesh 
enriched to quadratic elements. It can be observed that the error indicator correctly picked 
up the shock as the maximum error region. It is important to note that figure 4.2.6 presents 
effectivity index 7 *, not the error indicator. Due to the presence of the value of error in the 
denominator of the definition of 7 k (4.37), the effectivity index will often exhibit overshoots 
in the areas of low error (division by small numbers). That explains presence of high values 
of effectivity index in front of the low shock or in front of the plate in the next example. 
The global effectivity index for this problem was 7 = 7.7 and a standard deviation of local 

effectivity indices a = 1.67. 


Example 2: Viscous Flow Over a Flat Plate 

The two-step algorithm was also used to model the viscous flow past a flat plate. The 
problem being modeled was designed by the following data: 

• Mach number, Af = 3 

• Reynolds number, Re — 500 

• Free stream temperature ^ = 80° A' 

• The temperature of the plate, T w = 228°A 
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MIN=0.0012359 
MAX=0.9332073 
ERR0R=2. 14 18225 
D.O.F= 289 


Figure 4.2: (a) Flow over a blunt body. Distribution of error indicators. 



MIN=0. 166661 
MAX*8 .4579004 

D.O.F= 289 


Figure 4.2: (b) Flow over a blunt body. Local effectivity indices. 


92 



The finite element mesh is shown in Fig. 4.3. We applied initial h and p refinements to 
introduce appropriate layers of small higher order (up to p = 3) elements along the plate 
to resolve the boundary layer phenomena. Different shades of gray in Fig. 4.3 correspond 
to different orders of approximation. Elements with only their sides shaded are anisotropic 
elements with a higher order approximation in the direction perpendicular to the plate only. 
The solution of the flat plate problem in terms of density contours is presented in Fig. 4.4. 

Since the two-step algorithm consists of three linear steps, we performed an error estima- 
tion for each step. Similarly, as in Example 1, the exact errors involved in effectivity indices 
analysis were replaced by the errors obtained as differences between the actual solutions of 
Euler, momentum and energy steps, and the corresponding solutions obtained by enriching 
the order of approximation by 1 throughout the mesh, and performing one Euler or momen- 
tum, or energy time step, respectively. These differences were then measured in the energy 
norms defined by the bilinear forms associated with these steps, symmetrized as described 
in previous sections. 

Figures 4.5, 4.6, and 4.7 present distributions of the error indicators and local effectivity 
indices for the three steps of the two-step algorithm. The global effectivity indices 7 and 
standard deviations of local effectivity indices, 7F, in this problem were as follows: 


Euler step 

06 

II 

, ~G = 6.2 

momentum step 

7 = 25.9 

, a = 5.8 

energy step 

oo 

CO 

II 

II 

It* 


4.1.3 Relative Error Estimate 

The idea of the relative error estimate is to compare the finite element solution on a current 
mesh with a solution obtained on an enriched mesh and to measure the difference between 
the two solutions in a suitable norm. The enrichment of the mesh is done by raising the 
order of approximation of all elements by one. Of course, solving the problem on the enriched 
mesh is much more expensive than obtaining the original solution, so the method apparently 
does not seem very reasonable. However, if the original solution is a result of some expensive 
iterative process (such as, for instance, converging to a steady state solution in the case of 
viscous flow problems), then performing a single extra linear step on an enriched mesh is not 
a significant part of the total cost of the computations. In addition, solving of this problem 
can be performed with an iterative equation solver with a very good initial guess and with 
a very limited number of iterations (limited even to just one iteration). 

As a norm measuring the difference between the two solutions, one can use the energy 
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MIN= 0.389E-05 
MAX=0.0302722 
ERROR=0.0961987 
D.O.F= 958 


Figure 4.5: (a) Flat plate problem. Error indicators for the Euler step. 




MIN=0.0557633 
MAX=1 1.438621 

D.O.F= 396 


Figure 4.5: (b) Flat plate problem. 


Local effectivity indices for the Euler step. 
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MIN=0.0472422 
MAX= 12.23076 

D.O.F= 396 


Figure 4.6: (b) Flat plate problem. Local effectivity indices for the momentum step. 
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MIN=0. 12675 
MAX=29.698 

D.O.F= 396 


Figure 4.7: (b) Flat plate problem. Local effectivity indices for the energy step. 
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norm in the case of symmetric problems, or the norm defined by the symmetrized bilin- 
ear form of the original problem in case it is nonsymmetric but symmetrizable. With such 
choices of norms, the element residual method discussed in the previous section is an approx- 
imation of the relative error estimate. In fact, the element residual method approximates 
errors defined by the relative error estimate by expressing them in the form of local element 
contributions which are evaluated without actually solving the problem on an enriched mesh. 


4.2 Directional Adaptation Indicator 

The error indicators calculated from the element residual method have been used success- 
fully for the h- refinement for a number of hypersonic inviscid and viscous problems. Dur- 
ing the last year of this project, we have also implemented directionally-dependent error 
estimate schemes applicable to the h-p compressible flow solver. These directional adap- 
tation indicators will be discussed in this section. The current h-p data structure allows 
two kinds of mesh adaptation: h-refinement (refine/unrefine elements) and p-enrichment 
(isotropically /ansitropically increase the spectral orders of elements). Although the present 
h-p data structure only allows directional p-enrichment, the methodology discussed here is 
applicable to both directional ^-refinement and p-enrichment in two- and three-dimensional 
problems. 

It should be noted here that, in general, there exist no formal definition of directional error 
estimate - error norms used in the adaption process are defined in a full three-dimensional 
or two-dimensional spaces. The goal of our research is to provide directional adaptation 
indicator, which can choose an optimal refinement /enrichment direction. By optimal we 
understand a direction which provides maximum reduction of error norm due to a directional 
refinement / enrichment . 

According to this definition, the most natural way of defining directional adaptation 
indicator would consist of the following steps: 

1. try to refine/enrich the element in each of master directions (two or three depending 
on problem dimensionality), 

2. for each trial direction, estimate the error after the refinement, 

3. choose the adaptation direction which provides greatest error reduction. 

The above method, although formally correct, would be computationally too expensive. 
For practical purposes, we adopt two approaches which provide the directional adaptation 
indicator as a relatively simple and inexpensive extension of basic error estimation proce- 
dures. Construction of such an indicator is presented below. For the sake of clarity, we focus 
on a two-dimensional case. Extensions to three dimensions are immediate. 
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The first approach is based on the element residual method. Recall that the error in- 
dicator function <p K is computed based on the element enriched space, X h,p+\{K). This 
function can be effectively used as an indicator to determine the directionality of the error 
for that element. A natural choice is to use the norm of the directional derivative of the 
error indicator function in each coordinate on the master element 

\Mlx = J nK (^? dx * nd IkJliU = (4 ' 39) 

as the directional indicator. The actual refinement /enrichment direction is that of maximum 
norm of derivative of error function. This procedure is rather intuitive and theoretically 
unexplored, however, it has received a consistent support among researchers in the area 
of error estimation. The effectiveness of this directional adaptation indicator can only be 
confirmed by numerical experiment. 

The second approach is based more on interpolation error estimator. In particular, it 
focuses on different contributing components of the semi-norm of the solution: 


\U\ljc = \\Ut\\ljc + \\U, n \\l <K 

For practical purposes, the exact solution u can be replaced with the finite element 
solution U. Then, to determine the possible enrichment directions to improve the quality of 
the solution, one can utilize the norms of directional derivatives of the solution 

= L { ^ rjx {4M] 


Note that these values are only the local properties - they represent, for each element, the 
directional variations of either the error function or the solution. By selecting one of these 
norms and nnrm a.1iy. ing the ^-derivatives with respect to the sum of the two derivatives, a 
directional adaptation indicator can be defined as 




(4.41) 


or 



(4.42) 


The normalization gives 4>k v & value between 0 and 1, with the “indication” of enrichment 
in £- or 17 -direction according to whether the values is close to 0 or 1 , respectively. 


In practice, we first use the element residual method presented in the previous section to 
determine whether an element needs to be enriched. If the residual error for an element ex- 
ceeds the user specified threshold value, we then compute the direction adaptation indicator 
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defined above. According the user selected values for 61 and 62, where 0 < 61 < 62 < 1, the 
element is enriched as follows: 

• anisotropically enrich in ^-direction if 0 < <f>K V < f>i 

• anisotropically enrich in ^-direction if 62 < ^Kr, < 1 

• isotropically enrich in both directions if 61 < 4>Kr> < 62 

The values of 61 and 62, which control the directionality of the p-enrichment, are currently 
being selected based on numerical experience with typical values ranging between 0.2 through 
0.4 for 61 and 0.6 through 0.8 for b?. 

Numerical Example: Carter's Flat Plate Problem 

The directional adaptation indicator described above has been applied to the Carter’s flat 
plate problem. For this problem we have converged the solution on an initial linear graded 
mesh with two levels of 6-refinements as shown in Fig. 4.8. The corresponding density 
contours are also shown in Fig. 4.9. A map of the error indicator, 6k, calculated by the 
element residual method (as presented in Sec. 4.1) is shown in Fig. 4.10. Note, only the 
elements with error indicator values, 6k, greater than the threshold value (10 5 for this 
example) are considered for enrichment. Fig. 4.11 shows the corresponding plot of the 
directional adaptation indicator, 4>Kn, calculated from the directional derivatives of the error 
indicator function in equation (4.41). Note that the elements with 6k less than 10 5 (the 
ones not to be enriched) are not shaded in the plot. The next two figures, 4.12 and 4.13, 
present the resulting meshes after one p-enrichment pass for the values of (61,62) set to 
(0.2, 0.8) and (0.3, 0.7), respectively. 

Similar estimations were performed using the directional adaption indicators based on 
solution gradients. The corresponding plots of directionality indicator ipKn and the corre- 
sponding enriched meshes are shown in figures 4.14 to 4.16, respectively. 

A careful study of the above results clearly shows that both proposed directional indi- 
cators perform a good job of suggesting directional p-enrichment: isotropic at the plate tip 
region and normal to the wall within the developed boundary layer. It seems that directional 
indicators based on the interpolation error estimate (solution gradients) behaves in a more 
consistent fashion than directional indicator, based on residual error estimate. This is prob- 
ably because the solution of the local residual error estimation problem is more sensitive to 
element size and boundary conditions than the actual solution of the problem. 

Note that the well-known low effectivity of interpolation error estimators is not a problem 
in this case, because we are using the directionality indicator to choose between different 
enrichment directions, and not to estimate the actual directional error. 
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PROJECT: cartl_ie 


-MESH- 


PHLOW-C/2 



D.O.F=742 


Figure 4.8: Carter’s flat plate problem: Mesh with second level refinement. 
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PROJECT: cartl Je 


ERROR MAP 


PHLOW-C/2D 



enrichment threshold 

9 k = IQ' 5 


9k 


MIN=0.274E-06 
MAX=0.001255 
ERROR NORM=0.0 
D.O.F=742 


Figure 4.10: Map of the error indicators based on the element residual method. 
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PROJECT: cartl_ie 


DIRECTIONAL ADAPTATION INDICATOR 


PHLOW-C/2D 



6 k < 10- 5 
no enrich 


£ enrich 


£ and rj enrich j V enrich 


MIN=-0.1 
MAX =0.9123463 
D.O.F=742 




4>Kv 




Figure 4.11: Map of the directional adaptation indicator based on the residual error function. 
(User specified values of 6i and control the direction of p-enrichment) 
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PROJECT: cartl_ie 


-MESH- 


PHLOW-C/2D 



D.O.F=2258 


Figure 4.12: Enriched mesh according to residual-based <f>K v : 61 = 0 . 2,62 = 0.8 
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PROJECT: cartlje - MESH - 



1 2 3 4 5 6 7 8 


PHLOW-C/2D 


D.O.F=1977 


Figure 4.13: Enriched mesh according to residual-based 4>Kj)' = 0 . 3,62 = 0.7 
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PROJECT: caitl_ie 


DIRECTIONAL ADAPTATION INDICATOR 


PHLOW-C/ 



MIN=-0.1 

MAX=0.99891^ 

D.O.F=742 


Figure 4.14: Map of the directional adaptation indicator based on the solution gradients. 
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PHLOW-C/2D 



D.O.F=2080 


Figure 4.15: Enriched mesh according to fa,, based on solution gradient: bi = 0.2, 62 = 0.8 
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PROJECT: cartl_ie 


-MESH- 


PHLOW-C/2D 



D.O.F=1816 


Figure 4.16: Enriched mesh according to <f>Kn based on solution gradients: &i = 0.3, 62 = 0.7 
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4.3 Adaptive Strategies 

Once the distribution of errors is known, the decision must be made as to which elements 
should be refined. The h-p data structure allows two kinds of refinement: breaking elements 
(h-refinement) and enriching their spectral orders (p-refinement). As mentioned previously, 
the best improvement of the accuracy of finite element approximations can be achieved 
if a combination h- and p- refinements are performed. The problem of an optimal choice 
between h and p is, in general, still an open question. In the case of viscous flow problems, 
numerical experience suggests that p-enrichment is most desirable in boundary layer zones 
while in the rest of the computational domain the best choice is a linear approximation 
with h-refinements. The procedure of adapting meshes that we use in practice consists of 
performing exclusively h-refinements for steady state solutions until all the flow features have 
been resolved. Then, at the last stage of solving the problem, we enrich elements along the 
solid wall boundaries to provide high resolution of the viscous features of the flow, such as 
the skin friction and the heat flux. 

As an adaptive strategy for h-refinements we adopted the well known strategy of equili- 
brating errors which leads to optimal meshes in the case of elliptic problems. The procedure 
of adapting a mesh based on this strategy consists of the following steps: 


1. Integrate the solution forward in time until steady state is reached. 

2. Determine the distribution of the error indicators ej, i = 1, . . . , N, N being the number 
of elements. 

3. Refine the elements with the error larger than a certain percentage ai of the maximum 
error e mi : 

e,' > Cti • Cmu 

and unrefine the elements with errors less than ocj • 

4. Go to 1. 


The parameters a i and qj are user-specified and their values are assigned based on 
numerical experience. One practical selection of ai and a? is based on a percentage of the 
current number of elements to be unrefined and refined respectively. The general guideline 
is that in the case of flow problems, one should be rather generous with refinements and 
conservative with unrefinements. If smaller (refined) elements do not entirely cover shock 
regions, e.g. a shock across elements with different level of h-refinements, the small-large 
interelements may result in a misdisplacement of the shock. 
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5 Implicit /Explicit Procedures 


During the first extension phase of this project an adaptive implicit/explicit method for the 
solution of viscous compressible flow has been implemented. This method is being used 
mainly to increase the efficiency of the h-p adaptive Navier-Stokes solver. The algorithm is 
based on the general family of Taylor-Galerkin models, with several parameters controlling 
the actual implicitness of the scheme, and is combined with h-p mesh adaptation and adaptive 
selection of implicit and explicit zones. Several criteria for the selection of these zones have 
also been studied. The theoretical formulation of the implicit/explicit method has been 
accomplished previously (see reference [31]), and some numerical results have been obtained 
using a different code with ^-refinement capability. In the present work, the implementation 
of the implicit/explicit method is based on this experience, and the major task is to extend 
it for ftp-adaptive procedure. 

The basic idea of implicit/explicit algorithms is to combine the two methods to take 
advantage of the superior features of each. The major advantage of the explicit method is that 
element computations are relatively cheap and simple. Unfortunately this method suffers 
from stability limitations of the time step, which in some problems leads to prohibitively 
large numbers of time steps. 

The implicit algorithm allows for an application of larger time steps than the explicit 
method. Moreover, due to the existence of implicit boundary terms, it offers easy and 
straightforward control of natural boundary conditions, particularly those involving the vis- 
cous fluxes. An additional advantage is that with larger time steps no explicit artificial 
dissipation is necessary, which is very important in the calculation of boundary fluxes, par- 
ticularly wall heating rates. The major disadvantage of implicit methods is a much higher 
cost of element operations and a more complex and expensive solution of the resulting system 
of equations. 

In this section, the formulation and numerical implementation of an adaptive implicit/ 
explicit algorithm for compressible flows is presented. The algorithm is based on the general 
family of Taylor-Galerkin methods discussed in Section 2. 

5.1 Formulation of implicit /explicit schemes 

The algorithm of the implicit /explicit scheme must be designed so as to preserve stability, the 
conservative properties, and the required order of approximation. The procedure described 
below has been implemented in the code. 

We begin by partitioning the domain Q into subdomains f and SlW where explicit 
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and implicit schemes are to be applied, respectively, and 


n< E 'nn( ; ) = r £/ , ft (iJ > u ft (/) = ft 


It is convenient to assume that the interface between the two regions coincides with the 
element boundaries. 

It can easily be observed that the differential equations to be solved on the two subregions 
are different due to different implicitness parameters applied in each zone: 7 ^ in 

the implicit zone and a^ E \^ E \ 7 (E) = 0 in the explicit zone . Therefore, the variational 
formulation (2.41), based on the assumption of constant implicitness parameters, cannot 
be applied to the domain ft. Instead, it can be applied separately to each subdomain with 
additional continuation conditions across the interface. These conditions represent continuity 
of the solution and satisfaction of the conservation laws across the interface and are of the 
form: 


u (E) 

= u<») 

F (E)° 

= jf )c 

M E) 

= a!'» 

p(E)V 

= ri' ,v 


on Tei 


(5.1) 


where index n refers to the outward normal for the corresonding region (n* E ) = — n^). The 
continuity requirement also pertains to the test function, so that = V . Note 

that for general weak solutions of Euler equations the solution U need not be continuous 
across the interface. However, for regularized problems and finite element interpolation, the 
continuity of U is actually satisfied. 

If the variational statement is formulated for this problem, then in addition to interior 
integrals for each subdomain and regular boundary integrals, jump integrals across the in- 
terface appear in the formulation. These additional interface terms on the right-hand side 
are of the form: 



-At [F[ I)Cn + jrt £)Cn ] • V + At [Fi /)Vn + F {E)Vn ] • V 

+^! [(l - 2cP>) AWJ C ' + (l - 2a ,E >) A^F^"] . V* 


On the left-hand side of the variational formulation additional interface terms are of the 
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form: 


[ At [a<'U!' 1 ”Au''»n!'> + a< £ Ui E >"AL7 (E >nj £) l • V 
Jr IE 1 1 

- At [ 7 <'> (fi!'>AC/”»n!'» + Jf> ■ Al/<'>nf>) + T < £ > «>A u «/'n! E| + Pf >Al/< £ >nP>)] • V 

- ^ [(l - 2 q ( 'I) pWA'P A^AV^nf + (l - 2a ,E) ) ^* E U[ E) /l‘ £) At7lf ) n! E) ] ■ V* 

or, after reinterpretation of the linearized terms: 

/ At [ a <'>AF<„'> c + o< £ >AFi E > c ] V -At y»AF^ r - 7 < £ >AFi £ > v ] • V 
JT xe 

_ [(1 _ 2aW)pV) A yAF$ c + (1 - 2 aW)pWA[ E) AFW c \ . Vds 

In order to enforce interface conditions, the values of consecutive terms in the above formulas 
should be prescribed using equation (5.1). In this way, two first terms on the right hand 
side are set to zero. However, for other terms the direct application of interface conditions 
is impossible because of different coefficients for implicit and explicit components. In the 
present implementation of this scheme we set all the interface components to zero. This 
procedure preserves the continuity of fluxes and time accuracy across the interface up to the 
first order. Note that the artificial dissipation contributions, not presented here, are handled 
in exactly the same way as the natural viscosity. 

5.2 Selection of implicit and explicit zones 

The basic criterion for selection of implicit and explicit zones is simple: for a given time step 
all nodes which violate the stability criterion for an explicit scheme should be treated with 
the implicit scheme. According to this criterion, several options for an automatic adaptive 
selection of implicit /explicit zones were implemented: 

1. User- prescribed time step At: 

Within this option, the user prescribes the time step. All nodes satisfying stability 
criterion for the explicit scheme (with s certain safety factor) are explicit. This means 
that all the elements connected to these nodes are treated with the explicit scheme. 
On all other elements the implicit scheme is applied. 

2. Prescribed maximum CFL number: 

In this option, the user prescribes the maximum CFL number that can occur for 
elements in fl. The time step is automatically selected as the maximum step satisfying 
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this condition. The choice of a maximum CFL number may be suggested by the time 
accuracy arguments or the quality of results. 

3. A prescribed percentage of the domain is implicit. 

In this version, the user specifies the fraction of the domain which is to be treated im- 
plicitly. The elements with the strongest stability limitation (usually the smallest ones) 
are treated implicitly, the others are explicit. The time step is selected to guarantee 
stability of the explicit zone. 

4. Minimization of the cost of computations. 

This option has not been implemented in the code. It belongs to the advanced theory of 
the cost minimization, and it is included here to explicate some basic concepts which 
are crucial to the successful application of the implcit/ explicit procedure. In this option, 
the time step and the implicit/explicit subzones are selected to minimize the cost of 
advancing the solution in time (say one time unit). The algorithm is based on the fact 
that, for an increased time step, an increasing number of elements must be analyzed 
with the (expansive) implicit algorithm. The typical situation is presented in Fig. 5.1, 
which shows for different time steps the relative number of nodes that must be treated 
with the implicit scheme (to preserve stability). On the abscissa, the A tpE denotes the 
longest time step allowable for the fully explicit scheme (with certain safety factors). 
At fi denotes the shortest time step requiring a fully implicit procedure. The relative 
number of implicit nodes increases as a step function from zero for At < AtpE to one 
for At > tpi- Now assume that the ratio r of the computational cost of processing one 
implicit node to the cost of processing one explicit node is given. This ratio can be 
estimated relatively well by comparing the calculation time of element matrices and 
adding, for implicit nodes, a correction for the solution of the system of equations. Then 
the reduction of the cost of advancing the solution in time with the implicit/ explicit 
scheme, as compared to the fully explicit scheme, is given by the formula: 

*(A«) = ^(n< E ’ + rn<'>) 

Typical plots of the function R{At) are presented in Fig. 5.2. Shown here are the two 
cases: 

(a) the case of a small difference between fully explicit and fully implicit time steps— 
an almost uniform mesh 

(b) the case of a large difference between fully explicit and fully implicit time steps 
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Note that in either case, restrictions on the length of the time step should be applied, 
for example, from the maximum CFL condition. Otherwise the cheapest procedure 
would always be to reach the final time with one implicit step. 

From the plots in Fig. 5.2, the following observation can be made: for an essentially 
uniform mesh, the mixed implicit/explicit procedure does not provide savings of the 
computational cost — either a fully implicit or fully explicit scheme is the cheapest de- 
pending on the time step restriction. On the other hand, for very diverse mesh sizes 
the mixed procedure provides considerable savings. This means that the effectiveness 
of the mixed implicit/explicit scheme will be the best for large-scale computations with 
both very large and very small elements present in the domain. Some introductory nu- 
merical results confirming this observation are presented in Section 5.4. In the practical 
implementation of this method, the approximation of the function R(At) is automat- 
ically estimated for a given mesh. Then, the time step corresponding to the smallest 
R(At) is selected automatically (subject to additional constraints, in particular the 
CFLmax constraint). 

In our hp mesh adaptation procedure, all high order elements (p > 2) are handled by 
implicit scheme, i.e., all high order nodes are treated as implicit nodes regardless of the 
stability limitation. This is because, with hierarchial shape functions, and p < 2, even 
explicit procedure generates out-of-diagonal entries in the mass matrix. Since presently there 
are no efficient methods of lumping such a mass matrix, the explicit algorithm does not really 
improve efficiency for higher p orders. Therefore only the nodes of the linear elements are 
treated as candidates for implicit/explicit selection procedure. This implementation has the 
following two advantages: (1) since the high order elements are often used within boundary 
layers, application of implicit scheme in these regions provides faster convergence of the 
boundary fluxes and offers direct control of the natural boundary conditions; and (2) with 
the exclusion of high order elements from the explicit zone, the mass lumping procedure can 
be easily implemented and offers a very low computation cost in this region. 
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Figure 5.2: Reduction of the cost of computations due to implicit-explicit procedure. 
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The criteria described above axe based purely on a stability analysis. The calculation of 
the stable explicit time step for each linear node is based on the most conservative estimation: 
the minimum value is taken from all the elements connected to this node and also the elements 
whose constrained nodes are computed from this node. Nodes are defined as implicit or 
explicit according to whether the explicit stability criterion is violated or not. Once the 
implicit and explicit nodes have been defined, implicit and explicit elements are selected. All 
linear elements which contain at least one explicit node are defined as explicit, the remaining 
elements become implicit. 


5.3 Computational procedure 

In our compressible flow solver, the implicit /explicit procedure is combined with an h- 
p adaptation scheme. The range of implicit /explicit zones is redefined after every mesh 
adaptation and between adaptation after every prescribed number of steps. At every time 
step, the implicit /explicit procedure results in a system of equations 

(M + K (I) )U n+1 = R 

where the stiffness matrix has non-zero entries only for degrees of freedom in the implicit 
zone or for nodes with penalty-enforced boundary conditions. The solution procedure for 
the above system depends upon whether the consistent or lumped matrix is used, and the 
capability of the solver used for the system of equations, for example, whether it can solve 
efficiently a system with block diagonal structure. 

Based on the user specified implicit /explicit zone selection criterion described in previous 
section, the current version of the implicit/explicit solver uses a two-pass procedure: the 
first pass is a loop through explicit elements to solve for the solution of explicit nodes by the 
lumped mass method, and the second pass solves the remaining implicit degrees of freedom 
with a direct frontal solver. Incorporating this algorithm within the h-p data structure, the 
procedure can be listed as follows: 

FIRST PASS 

1 Loop through the explicit elements: 

1.1 form the unconstrained consistent mass matrix and right-hand-side (RHS) 
vector 

1.2 modify the mass matrix and RHS vector for contrained nodes 

1.3 lump the mass matrix and assemble RHS vector for explicit nodes 

2 Obtain the solution for explicit nodes 
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SECOND PASS 


3 Loop through all elements: 

3.1 For explicit elements containing implicit nodes: 

3.1.1 form the unconstrained mass matrix and RHS vector 

3.1.2 modify the unconstrained mass matrix and RHS vector for constrained 
nodes 

3.1.3 lump the mass matrix 

3.1.4 apply boundary conditions 

3.1.5 for implicit nodes, store the stiffness matrix and RHS vector in a frontal 
solver format 

3.2 For implicit elements: 

Following the standard procedure, form the stiffness matrix and RHS vector, 

and then store them in frontal solver format 

4 Obtain the solution for the implicit degrees of freedom using direct solver. 

Remarks 

1. The calculation of mass matrix in step (1.1) needs to be done only for a single solution 
component. 

2. If the global lumped mass matrix is saved, it needs to be recomputed only after the 
mesh is adapted or the implicit /explicit zones are redefined. 

3. The lumped mass matrix is obtained by summing the rows of the “constrained” consis- 
tent mass matrix. For boundary nodes in explixit elements, the lumping is performed 
before the application of the boundary conditions. 

4. If the application of penalty-enforced boundary conditions produce off-diagonal entries 
to the stiffness matrix (e.g. solid- wall and no-flow boundary conditions), the boundary 
nodes are treated as implicit nodes and their solutions are solved in the second pass. 

Numerical Illustration of Implicit/Explicit Procedures 

The first example is the Carter flat plate problem solved on a mesh with 3 levels of h- 
adaptation. The density contours are shown in Fig. 5.3. Figs. 5.4a and 5.4b represent the 
meshes when 40 % and 80 % of the degrees of freedom are selected as implicit, respectively. 
The implicit elements are shaded dark, the explicit elements containing at least one implicit 
node or node on the boundary which gives off-diagonal entries are shaded light, and the 
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explicit elements containing only explicit nodes are not shaded. A complete set of timing test 
runs has also been performed on this problem. These test used the third implicit/explicit zone 
selection option — the user specifies percentage of the domain to be treated implicitly. The 
results are tabulated in Table 5.1, showing the percentage of implicit zone, the timestep size, 
At, the equivalent CFL number (based on fully explicit sheme with safety factor set to 1). 
and the cost reduction factor (also based on fully explicit scheme). When the implicit/explicit 
zone is selected based on the cost minimization with the limitation on the time step controlled 
by the condition CFL < 10, about 65 % of the nodes are treated as implicit nodes. The 
computation cost of reaching the steady-state solution is reduced by a factor of 2 with respect 
to the fully explicit algorithm. 

The second example is a blunt body with incident shock (the detailed description of the 
problem is presented in Section 7). A relative coarse initial mesh consisting of 18 by 40 
linear elements is used for the inviscid solution. The corresponding pressure contours are 
shown in Fig. 5.5. The stable time steps of each node for the explicit scheme (based on 
CLF = 1) are contoured in Fig. 5.6 and range from 0.0326 to 0.3408. The implicit and 
explicit zones selected for a given time step of 0.1 are plotted in Fig. 5.7. The corresponding 
cost reduction factor is about 0.71. The same mesh with 2-levels of ^-adaptation for the 
viscous solution and the pressure contour are shown in Figs. 5.8 and 5.9, respectively. The 
stable time steps of the nodes for explicit scheme are in the range between 0.006727 and 
0.3412. The implicit and explicit zones selected, shown in Fig. 5.10, are based on the second 
option with a specified maximum CFL of 3. For visual clarity, the drawing of element grid 
lines is suppressed and the colormap is adjusted so that the light, medium and dark shades 
represent the fully explicit elements, explicit elements containing implicit nodes, and fully 
implicit elements, respectively. The nodes treated implicitly cover about 37 percent of the 
domain, and cost reduction factor with respect to the fully explicit method is about 0.92. 
Compared to a fully implicit method (which was previously the only method available in this 
project) the cost reduction factor of the implicit/explicit method is about 0.14. It should 
be noted that, although the pressure contours show reasonable resolution of shocks, the 
clustering of the elements near the blunt body is still too coarse to resolve the flow features 
in the viscous layer. The numerical results and discussion using a finer mesh are presented 
in Section 7. 
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PHLOW-C/2D 
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Figure 5.4: Implicit and explicit zones for the Carter fiat plate problem (a) with 40 % of 
nodes treated implicitly and CFL S* 4.5 and (b) with 80 % of nodes treated implicitly and 


CFL S 29. 


125 



%of 

Implicit 

At 

Equivalent 
CFL number 

Cost 

Reduction 

0 

.000335 (A t FE ) 

1 

1 

10 

.000446 

1.33 

0.97 

20 

.000755 

2.25 

0.73 

30 

.000906 

2.87 

0.81 

40 

.001511 

4.51 

0.63 

50 

.002012 

6.01 

0.62 

60 

.002927 

8.74 

0.55 

70 

.004415 

13.18 

0.45 

80 

.009791 

29.23 

0.25 

90 

.011064 

33.02 

0.27 

100 

.1124 (A t F! ) 

(335.5) 



Table 5.1 
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Figure 5.6: Contours of stable time-step sizes of nodes. 
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Figure 5.7: Implicit and explicit zones selected by specifying At = 0.1 
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PROJECT: bb4_r 


-MESH- 
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Figure 5.8: Two level fc-adapted mesh for viscous solution. 


D.O.F=4989 
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PRESSURE 


PHLOW-C/2D 
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MAX=6.6864227 


Figure 5.9: Pressure contours of viscous solution. 
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max CFL = 3. 



6 Some Nonstandard Algorithms 


In this section we present several nonstandard finite element procedures which are necessary 
to efficiently implement the h-p finite element methodology. 


6.1 Integration and Underintegration Procedures 

Solving the boundary value problems associated with the Taylor- Galerkin approximation 
and the approximation of the viscous step involves enforcing essential and mixed boundary 
conditions, such as the no-penetration boundary condition on solid walls. These boundary 
conditions are currently being enforced by means of a penalty function. It is well known 
that this approach requires special integration procedures guaranteeing that the number of 
integration points over any patch of elements on the boundary matches the number of degrees 
of freedom in this region. Otherwise, if the number of integration points exceeds the number 
of degrees of freedom, a locking effect may result when enforcing more boundary constraints 
than needed. Practically, this means that integration points must be associated with nodal 
locations and their number must be equal to the number of shape functions corresponding 
to a given node. 

In two dimensions the usual way of dealing with this problem is to introduce appropriate 
underintegration on the boundary with mixed conditions, i.e., integration which for a p-th 
order one-dimensional boundary of the element introduces p + 1 points with two of them 
at the endpoints. The Gauss-Lobbato integration scheme satisfies such conditions and is 
frequently used. 

In three-dimensional h-p finite element methods, the situation is more complicated, since 
one may now encounter constrained nodes and nonuniform distributions of degrees of ap- 
proximation p on the two-dimensional surface of the computational domian. Consider, for 
example, the situation presented in Fig. 6.1, where both active and constrained nodes are 
shown. If we take the element K \ , unconstrained and with uniform orders p, then a tensor 
product of p + 1 one-dimensional Gauss- Lobatto integrations will have the proper number 
and distribution of integration points. If we consider, however, an element with nonuniform 
p’s, like K 2 , then no tensor product of one-dimensional procedures will be suitable as such 
products distribute integration points according to rectangular patterns not adequate for 
nonuniform distributions of degrees of freedom. The situation is more complicated for con- 
strained elements, like K 3 : integration points along the side adjacent to element K+ must 
be chosen such that they coincide with integration points of the larger unconstrained neigh- 
bor. Otherwise, more conditions on the global shape functions would be imposed over both 
elements than the number of degrees of freedom. Figure 6.2 presents a possible distribution 
of integration points which would be satisfactory from the point of view of avoiding locking 
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effects. 

The problem that emerges is to construct integration procedures which are accurate 
enough but which also satisfy the indicated restrictions on the numbers and locations of 
sam pling points. To avoid this problem, note that the integration of penalty terms does not 
require any accuracy of the integration procedure: the penalty approach is an implicit way 
of imposing pointwise conditions on the solution (collocation). The only reason for using 
some specific integration procedures (i.e., with some accuracy) is that not only are penalty 
terms integrated on the boundary, but also other quantities (like stresses, etc.) are computed 
on the boundary which we wish to integrate with sufficient accuracy as well. This suggests 
that we simply use separate procedures for integrating penalty terms and all the remaining 
quantities. In the first case, use sampling points (not necessarily integration points) whose 
number and location corresponds to locations of nodes and their orders p, with (say) unit 
wei gh ts. In the second case, use ordinary numerical integration procedures, for instance 
Gaussian quadratures. An example of a choice of sampling points for integrating penalty 
terms is shown in Fig. 6.3: they are distributed uniformly along unconstrained sides and 
the interior of the element wall with their number being p + 1 or (p + l) 2 , respectively. For 
constrained sides, sampling points are chosen to coincide with those of the larger neighbor 
for which the common side is not constrained. 
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Figure 6.2: Integration point locations used to avoid locking. 
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The problem of enforcing mixed boundary conditions involves also a choice of approxima- 
tion of the outward normal to the boundary. The simplest solution is to define the outward 
normal at a point common to several elements as an average of outward normals evaluated 
for all these elements. 


6.2 Routines for Performing Refinements in Boundary Layer 
Zone 

As mentioned in the section on adaptivity, p-refinements along the solid wall boundaries are 
currently not performed automatically, but rather “by hand.” That is, a user decides which 
elements should be enriched. In fact, designing a mesh in the boundary layer zone may also 
require some other operations like breaking or unrefining some elements. If one solves a test 
problem with a relatively small number of elements, such an operation can be performed 
interactively. In the case of a real life problem, however, the user would have to indicate a 
large number of elements to be refined or unrefined. 

A method has been developed that significantly reduces the effort of interactive generation 
(refinement and unrefinement) of simple meshes. The method consists of mapping a possibly 
curvilinear portion of a mesh into a rectangular domain with element sides parallel to the 
sides of a rectangle. Of course, the method can be applied only to pieces of the mesh for 
which the pattern of initial elements is topologically equivalent to a uniform rectangular 
mesh. 

The idea of using the mapped image of the mesh to perform the required refinements is 
simple. First one identifies the rectangular coordinates of a group of elements that they wish 
to refine (this is an easy operation as mapped elements constitute rectangular patterns). 
Then one performs a specified refinement by giving the coordinates of rectangles covering 
the considered subdomain. For instance, one may prescribe: break elements between x = 0 
and x = 1 and between y = 0.5 and y — .75, etc. This way of generating meshes and 
simple refinements is still interactive yet it allows one to perform massive refinements with 
a minimal effort. (See the User Manual for additional details.) 

6.3 Postprocessing in Three Dimensions 

When solving three-dimensional problems, special attention must be devoted to the problem 
of displaying the solution. The three-dimensional finite element code is equipped with a 
postprocessing package which can display contour maps of the solution on the boundary and 
on the desired cross sections of the computational domain. This postprocessing is rather 
expensive and it is not always sufficient for displaying interesting features of the solution. 
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Figure 6.3: Sampling points foT integrating penalty terms. 
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For instance, contour maps give only a rather inaccurate image of quantities defined on 
the boundary of the domain such as, say, heat flux or skin friction. Three-dimensional 
perspectives of such functions give a much better idea of their behavior. 

Following this motivation an algorithm was developed for displaying three-dimensional 
perspectives of the solution dependent functions on the portions of the boundary of a three- 
dimensional domain. In the first step, the indicated portion of the boundary, possibly curvi- 
linear, is continuously mapped into a flat plane. Two kinds of mappings are possible; a 
projection in a given direction or the mapping “stretching” a two-dimensional mesh of dis- 
torted rectangles on the boundary into a plain mesh of distorted rectangles on the boundary 
into a plain mesh of rectangles, the transformation described in the previous section. In 
the second step, a three-dimensional perspective of the solution is displayed over the plane 
domain. 

The procedure is inexpensive if compared to contouring in three dimensions and it is very 
helpful in showing intricate features of a three-dimensional solution. 

6.4 Solution Correcting Procedures 

When integrating the Euler or Navier-Stokes equations one frequently encounters a cum- 
berson phenomenon: the solution evolves to a physically unacceptable state with negative 
densities or pressures. Such situations make further integration impossible unless one artifi- 
cially corrects the solution by “pushing it back” to physically meaningful values. 

A procedure for performing this operation was designed. It can be outlined as follows: 
First, one verifies if 


’ p > 0 

. * = e - ^ (m 2 + ml) fp > 0 

where p is density, t internal energy, m, n momentum components, and e is the total energy. 
Note that the above conditions imply that e > 0 and the pressure p = (7 — l)t > 0 . If either 
p < 0 or 1 < 0, one replaces p and e by new values p* and e* defined as a projection of the 
point (p, e) € R 2 onto a set A C £ 2 of physically acceptable densities and energies: 

A = j(p, e) € £ 2 : p > 0 and pe > i (m 2 + n 2 ) j 

The correcting procedure described here has proven very useful especially when starting the 
integration process when the initial solution may not satisfy given boundary conditions and 
lead immediately to nonphysical states. 
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7 Numerical Examples 


Several problems were solved to test the algorithms presented in the previous sections. The 
problems range from the simplest two-dimensional inviscid flow over a wedge to three- 
dimensional viscous flows with high Mach and Reynolds numbers. In all of the numerical 
examples presented here we have focused on providing solutions exclusively for steady state 
problems. The integration in time process is a method of obtaining such solutions. After 
converging to steady state, error estimation and mesh adaptation are performed according 
to rules presented in Section 4. The residual error indicators are used as the basis for mesh 
refinement with the equidistribution of errors as a refinement strategy. 


The stable time step 
formula 


A t used in the calculations has been evaluated by the following 


At = CFL min 


h_K 

M + c 


(7.1) 


where |«| is the maximum magnitude of the fluid velocity in the element K, c is the speed 
of sound, and CFL is the Courant, Friedrichs, Lewy stability coefficient [6]. In addition, all 
three types of artificial dissipation, the Lapidus and Morgan’s viscosities and a generalization 
of the last one for distorted elements have been used for various problems. The implicitness 
parameters and 7 in all the examples are set equal to one. As a linear equation solver we 
have used the block Jacobi algorithm accelerted with GMRES. 


Example 1: Inviscid Flow Over a Wedge on h- Adapted Meshes 

We consider the problem of supersonic flow over a wedge with the following prescribed data: 

Mach number M = 3.0 
Inclination of the wedge 6 = 20° 

The initial meshes consist of 12 x 5 linear, quadratic, and cubic elements. The Lapidus 
viscosity constants for these three cases were assumed as c* = 1, c* = 0.15, c* = 0.07, 
respectively. 

The ^-adaptive meshes are shown in Figs. 7.1, 7.3, and 7.5, and the corresponding density 
contours are given in Figs. 7.2, 7.4, and 7.6. Comparing these figures one finds that all three 
meshes have captured the shock within two or three elements and all of the shock angles are 
virtually identical. In addition, the density contours obtained from the linear mesh (Fig. 7.2) 
show a much tighter shock pattern than the quadratic or cubic meshes (Figs. 7.4 and 7.6). 
This is most likely a result of the two additional levels of ^-refinement and an introduction 
of approximately 50 percent more degrees of freedom in the shock region than in either of 
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the other two cases. Finally, note that all these cases show very little of any reflection of the 
shock at the subsonic outflow boundary on the upper surface. 
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D.O.F= 1507 


Figure 7.1: Flow over a wedge problem on an h-adaptive mesh of linear elements. Final 
mesh. 
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Figure 7.2: Flow over a wedge problem on an ^-adaptive mesh of linear elements. Density 
contours. 
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Example 2: Inviscid Flow Over a Blunt Body on h- Adapted Meshes 


As a second test problem, we considered the supersonic flow over a cylinder with Mach 
number M = 6.0. The problem was solved starting with initial meshes of 16 x 16 linear 
elements and 8x8 quadratic and cubic elements. The viscosity constants c* were set to 1.0, 
0.07, and 0.035, respectively. The final /i- adapted meshes are presented in Figs. 7.7, 7.9, and 
7.11, the corresponding density contours are shown in Figs. 7.8, 7.10, and 7.12. Comparing 
these results one finds again that all three meshes have predicted virtually identical stand-off 
distances for the bow shock, and have captured the shock within two or three elements. The 
resolution of the shock for the linear and quadratic meshes in this case is again somewhat 
better than for the cubic mesh, which is most likely due to the number of degrees of freedom 
introduced in the shock region. Finally, note that the density contours for the cubic elements 
are the smoothest of the three meshes. This same pattern can be observed in Figs. 7.2, 7.4, 
and 7.6 of Example 1. 
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Figure 7.7: Blunt body problem on an h-adaptive mesh of linear elements. Mesh after two 
levels of refinements. 
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Figure 7.10: Blunt body problem on an h-adaptive mesh of quadratic elements. Density 
contours. 
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Example 3: Carter’s Flat Plate Problem 

The third test case we considered was Carter’s Flat Plate Problem [5] with the following 
data 

Moo = 3, Re L = 1000, Pr = 0.72, 7 = 1-4, Too = 390 [°R] (7.2) 

The geometry and boundary conditions are shown in Fig. 7.13. The following flow featues 
are recognized: 

• a singularity exists near the leading edge of the plate, 

• a curved bow shock is developed from the tip of the plate, and 

• a boundary layer is formed along the plate. 

The flags in the figure (KIND = 1,. . . , 5) correspond to various boundary conditions incor- 
porated in the code. The solid wall temperature is 

T wa u = 1092[°it] (7.3) 


A nonuniform initial h-p finite element mesh, as shown in Fig. 7.14 with second and 
third order elements defined along the solid wall boundary, has been used to march to the 
steady state condition. In order to reduce the numbers of degrees of freedom in most of the 
boundary layer region (except near the stagnation point), anisotropic elements were used 
which have a high order of approximation in the direction normal to the plate and a linear 
approximation along the plate. 

Using the pure /i-adaptive strategy described in Section 4, the mesh shown in Fig. 7.15 
was obtained. The corresponding density contours are presented in Fig. 7.16. The next 
two sequences of adaptive meshes are shown in Figs. 7.17 and 7.19 with the corresponding 
density contours in Figs. 7.18 and 7.20. After the third adaptive refinement, the solution 
process was stopped when the error reached the prescribed error tolerance. In all steps, the 
CFL constant was set to 0.5. 

Figures 7.21 and 7.22 compare the computed profiles of pressure and skin friction dis- 
tributions along the solid wall with the original results of Carter, showing good agreement 
away from the stagnation point and the ability of the method to also resolve the stagnation 
point itself. Figure 7.23 shows the computed heat flux along the solid wall (not presented 
in Carter’s report). Note the positive sign in the vicinity of the tip of the plate (the gas 
heating plate). The stabilizing effect of higher order elements near the solid wall can also be 
observed in the pressure contours shown in Fig. 7.24 corresponding to the final mesh. 
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Figure 7.13: Carter’s flat plate problem. Geometry and boundary conditions. 
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Figure 7.14: Carter’s flat plate problem. A nommiform initial h-p mesh. 
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Figure 7.15: Carter’s flat plate problem. An optimal mesh with maximum level of refinement 
equal 2 (Mesh 1). 
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Figure 7.16: Carter’s problem. Mesh 1. Density contours. 
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Figure 7.17: Carter’s problem. Mesh 2. 
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Figure 7.18: Carter’s problem. Mesh 2. Density contours. 
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Figure 7.19: Carter’s problem. Mesh 3. 
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Figure 7.20: Carter’s problem. Mesh 3. Density contours. 
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Figure 7.21: Carter’s problem. Pressure profiles along the plate. 
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Example 4- Flat Plate Problem with Re = 10,000 

As a continuation of Example 3, we have also performed an analysis of the flat plate problem 
with a Reynolds number of one order of magnitude larger. This required solving example 
problem 3 in a domain roughly ten times larger than before. In order to avoid an excessive 
number of degrees of freedom in the initial meshes, resulting from higher order elements in 
the initial stage of the solution procedure, higher order elements were introduced only in the 
last adaptive step after all A-refinements were performed. 

Starting with an initial mesh of only 4x8 elements, a sequence of consecutive optimal 
linear meshes were obtained using the adaptive strategy discussed in Section 4. A total of 7 
refinements were generated, each computed only when the optimal mesh for a particular level 
of refinement had been obtained. The final mesh and the corresponding density contours 
are shown in Fig. 7.25. Figure 7.26 presents a three-dimensional perspective of the same 
density function showing a clear separation of the shock from the boundary layer. 

For coarse meshes (low levels of refinement), the computed boundary layer is primarily 
due to numerical viscosity, the corresponding viscous quantities, especially near the leading 
edge, being far from those obtained in the previous example. Figure 7.27 presents, for 
instance, the profile of heat flux along the solid wall (compare with Fig 7.23). 

The situation drastically changes, however, when a layer of second, third, and fourth 
order elements are introduced into the mesh (see Fig. 7.28). The corresponding heat flux 
profile, shown in Fig. 7.29, agrees qualitatively with the results obtained in the previous 
example. 
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Figure 7.25: Flat plate problem with Re = 10, 000. Density contours and optimal mesh of 
linear elements with the maximum level of element = 7. 
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Figure 7.26: Flat plate problem with Re = 10,000. A three-dimensional perspective of the 
density function on the optimal mesh of linear elements. 
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Figure 7.27: Flat plate problem with Re = 10,000. Heat flux coefficient profile along the 
plate for the mesh of linear elements. 
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Figure 7.28: Flat plate problem with Re = 10,000. Blowup of the optimal mesh enriched in 
the boundary layer. 
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Figure 7.29: Flat plate problem with Re = 10,000. Heat flux coefficient profile along the 
plate for the enriched mesh. 
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Example 5: The Holden Compression Comer Problem, M = 5, Re = 30, 000 

The next two-dimensional test case modeled was the viscous flow over a compression corner. 
The problem was solved for the following data: 

M =5 
Re = 30,000 
T wa ii = 288° K 
^ = 80° K 

angle of inclination 15° 

The initial mesh consisted of 7 x 17 linear elements. Three levels of h-refinements were 
subsequently performed leading to a good resolution of the shocks inside the computational 
domain. Finally, layers of quadratic, cubic, and fourth order anisotropic elements were in- 
troduced along the solid wall boundary. The elements in the neighborhood of the stagnation 
point were not enriched as we experienced some stability problems with higher order elements 
in this area. The final h-p mesh is shown in Fig. 7.30. 

The solution, contours of density and profiles of the density, velocity, temperature and 
pressure coefficient are shown in Figs. 7.31-7.43. As in the previous example a stabilizing 
effect of higher order approximations on the behavior or pressure can be observed. The 
profiles of the heat flux and the skin friction are shown in Figs. 7.44 and 7.45. A dramatic 
change in the quality of these two fluxes is observed as we move from the section of the 
boundary with linear elements only to the enriched part. 

The higher order approximation resulted also in a satisfactory resolution of another fine 
feature of the solution: the recirculation of the flow. This is shown in Figs. 7.46 and 7.47 
where the directions of the flow and the negative contours of Uj-velocity in the zone of 
recirculation are presented. 
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Figure 7.30: Holden’s compression corner problem, h-p adaptive mesh. 
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Figure 7.31: Holden’s compression corner problem. Density contours. 
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Figure 7.32: Holden’s compression corner problem. Profile of density along vert 1. 
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Figure 7.33: Holden’s compression corner problem. Profile of velocity along vert 1. 
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Figure 7.34: Holden’s compression corner problem. Profile of temperature along vert 1. 
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Figure 7.36: Holden’s compression corner problem. Profile of density along vert 2. 
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Figure 7.39: Holden’s compression corner problem. Profile of pressure coefficient along vert 
2 . 
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Figure 7.40: Holden’s compression corner problem. Profile of density along vert 3. 
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Figure 7.41: Holden’s compression corner problem. Profile of velocity along vert 3. 
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Figure 7.42: Holden’s compression corner problem. Profile of temperature along vert 3. 
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Figure 7.43: Holden’s compression corner problem. Profile of pressure coefficient along vert 
3. 
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Figure 7.44: Holden’s compression corner problem. Profile of heat flux coefficient along the 
plate. 
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Figure 7.45: Holden’s compression corner problem. Profile of skin friction coefficient along 
the plate. 
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Figure 7.46: Holden’s compression corner problem. Recirculation: velocity vectors near the 
corner. 
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Figure 7.47: Holden’s compression corner problem. Recirculation: contours of the uj com- 
ponent of the velocity. 
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Example 6: Inviscid Flow Past a 20° Wedge, M = 3 

Our first three-dimensional test case was the supersonic inviscid flow (M = 3) over a two- 
dimensional wedge (analyzed using the three-dimensional code). The initial mesh consisted 
of one layer of 4 x 8 linear elements. On the side surfaces we imposed symmetry boundary 
conditions, which enforce (weakly) the condition dU/dn = 0, which was intended to enforce 
the two-dimensional character of the flow. The solution (contours of density) and A-adapted 
mesh are shown in Fig. 7.48. One can observe that indeed the solution is essentially inde- 
pendent of the y-coordinate direction. Comparing this result with the two-dimensional case 
(Fig. 7.2), one observes a good qualitative agreement of the flow features and shock angle. 

Example 7: Inviscid Flow Around a Sphere, M = 6 

Next, the supersonic ( M = 6) flow over a spherical blunt body was solved. In the discretiza- 
tion, symmetry was enforced so that only one quarter of a half sphere is meshed. The initial 
mesh is generated by appropriately mapping a regular 7x7x4 rectangular mesh onto the 
quarter sphere. Such a mapping, if performed to cover exactly the octant of a sphere, can 
lead to severe distortions of some elements as it has to transform a rectangular domain into 
a topologically triangular domain; hence, the octant is not covered exactly, leaving some 
missing sections on the outflow boundary (see Fig. 7.49). On the lower and right planes of 
Fig. 7.49, symmetry boundary conditions are imposed. Modified Lapidus viscosity is used 
as the artificial dissipation mechanism in this problem. The h-adapted mesh and computed 
density contours axe shown in Figs. 7.50 and 7.51. The flow is characterized by a bow surface 
of the shock surrounding the body. 

Example 8: Viscous Flow Past a Flat Plate, M = 3 

The classical two-dimensional flat plate problem was also modeled using the three- 
dimensional code. The data for the problem are: M = 3, Re = 500, Too = 80 K, 
T wa u = 280° K. The computational domain was discretized initially by one layer of 4 x 8 
elements. After converging to a steady state solution the elements along the plate were 
refined and enriched to second and third order. Anisotropic p-enrichments are used, intro- 
ducing higher order approximations only in the direction perpendicular to the wall, so as 
to significantly improve the approximation of the boundary layer. Locally, the structure 
of the boundary layer is close to that of the one-dimensional case, the largest gradients 
being in a direction perpendicular to the wall. The enriched mesh is shown in Fig. 7.52, 
where anisotropically enriched elements are marked by shading their edges in the direction 
of enrichment. 
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Figure 7.48: Inviscid flow past a 20° wedge, M = 3. Density contours and an h - adapted 
mesh. 
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Figure 7.49: Flow around a sphere, M = 6. An initial mesh. 
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Figure 7.50: Flow around a sphere, M = 6. A-adaptive mesh. 
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Figure 7.51: Flow around a sphere, M = 6. Density contours. 
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After converging to a steady state solution on this mesh /i-adaptations were performed 
based on interpolation error indicators. The new mesh and density contours are shown in 
Figs. 7.53 and 7.54. Observe that major features of the flow, the shock and boundary layer 
zone, are well-resolved. Figs. 7.55 and 7.56 present plots of the skin friction coefficient and 
the heat flux coefficient. 
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Figure 7.52: Viscous flow over a flat plate, M = 3. h-p adapted mesh after one refinement 
iteration. 
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Figure 7.53: Viscous flow over a flat plate, M - 3. h-p adapted mesh after two refinement 
iterations. 
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Figure 7.54: Flat plate problem, M — 3. Re = 500. Density contours. 
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Figure 7.55: Carter flat plate, M = 3. Profile of skin friction coefficient along the plate. 




Figure 7.56: Carter flat plate, M = 3. Profile of the heat flux coefficient along the plate. 
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Benchmark Problems 

Example 9: Blunt Body Problem With a Type IV (Edney) Interaction 

The first benchmark problem is a supersonic viscous flow around a blunt body problem with 
an incident shock, as defined in Fig. 7.57 and described in [32]. During the first phase of 
this analysis we applied the adaptive finite element method to the inviscid case. A set of 
preliminary results for this case is presented here in Fig. 7.58 to Fig. 7.62. The initial mesh 
consists of 32 by 16 linear elements and is shown in Fig. 7.58. The final A-adapted mesh 
with 3 levels of refinement and the contours of density are shown in Fig. 7.59 and 7.60, 
respectively. The major deficiency in this solution is a poor resolution of the shock near 
the cylinder, where the error estimator (the residual technique) was not indicating large 
errors, at least as compared to the other shocks. To improve the quality of the solution, an 
additional manual mesh refinement was introduced in this region simultaneously with the 
automatic mesh adaptation procedure. A solution on this final mesh, Fig. 7.61, is shown in 
Fig. 7.62. (The above results was obtained by the two-step procedure described in Sec. 2.) 

During the second phase of the analysis, the work has focused on the viscous solution 
of the blunt body problem, where the freestream Reynolds number (based on the cylinder 
radius) was 2.00977* 10 5 , and the wall temperature was fixed at 530° R. Two different meshes 
have been used to solve this case. The first mesh (referred as mesh A), shown in Fig. 5.7, 
consists of 40 by 18 linear elements initially, and is graded by a geometric progression of 
factor 1.01 in the radial direction. To initialize the flow, the inviscid solution was calculated 
on a one- level uniformly refined mesh. This solution was used as a starting point for a viscous 
analysis. The mesh with two levels of A-adaptation and the pressure contours of the viscous 
solution are presented in Figs. 5.8 and 5.9. A 3D perspective view of the pressure is also 
shown in Fig. 7.63. While the result shows a reasonable shock resolution, the viscous features 
near the cylinder are still far from fully resolved and do not compare well with experiment 
data provided by NASA LaRC. This indicates a need to introducing more degrees of freedom 
along the cylinder to resolve the viscous phenomena. Notice that the smallest thickness of 
the element along the cylinder after two levels of refinement is still greater than 1% of the 
radius. Ref [32] also indicates that the dominant viscous region is less than 1% of the radial 
distance from the cylinder. 

In order to resolve viscous phenomena near the wall, a higher order approximation was 
introduced in the boundary layer. A mesh (referred to as mesh B) with such features is shown 
in Fig. 7.64. The boundary layer zone of width approximately 0.01 was covered by nine layers 
of second order elements with sizes geometrically graded toward the wall (with ratio q = 0.5). 
The size of the smallest element is was based on laminar boundary layer theory as well as on 
results presented in literature [32]. The rest of the computational domain is discretized by 
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40 by 18 linear elements. A final adapted mesh with two levels of ^-adaptation (with 7082 
degrees of freedom) is show in Fig. 7.65. The corresponding pressure contour and sonic line 
are presented in Figs. 7.66 and 7.67, respectively. The 3D perspective views of the density 
and pressure in the shock-boundary layer interaction area are shown in Figs. 7.68 and 7.69, 
where the embedded shocks in the jet zone can be clearly identified. 

The mesh with quadratic elements in the near-wall region for this case provided a much 
better resolution of viscous phenomena and favorable comparisons with experimental results. 
In particular, the comparison of pressure and heat flux along the cylinder with the exper- 
imental data are profiled in Figs. 7.70 and 7.71, respectively. The predicted wall pressure 
distribution is in very good agreement with the experiment results, except the small dif- 
ference of the shock impingement location. (The same discrepancy was also reported and 
explained in Ref [18], where it was essentially attributed to fluctuation in the experimen- 
tal flow pattern.) The plot of the heat flux shows a good correlation with experimental 
observations — better than many other numerical references with much finer meshes. Note 
that the oscillatory character of our numerical heat flux is a consequence of a temporary 
imperfection in our post-processing package — gradients are calculated and plotted at nodes, 
where their accuracy is the lowest. A better postprocessing algorithm would be to perform 
the gradient calculations at integration points and then project the solution to the nodes. 
Also note that a discrepancy in the magnitude of the heat flux distribution between the 
numerical and experimental data w as also reported in Refs. [18,32]. It should be mentioned 
that our computational grid has a much coarser circumferential resolution than those used in 
other references. The cylinder was covered by 40 quadratic elements in the circumferential 
direction (on the adapted mesh). Because of the curvature effects around the cylinder, high 
accuracy on such a coarse mesh can only be achieved by using high order elements. For 
example, mesh A, which consists of linear elements only, and with even more than 10,000 
degrees of freedom, failed to provide a reasonable resolution of heat flux. 

The numeri cal results on the last adapted mesh B were obtained using the im- 
plicit/explicit procedure. Because the stable time step sizes are always restricted by the 
extremely thin elements in the viscous layer, both fully implicit and fully explicit schemes 
are more expensive than the mixed implicit/explicit procedure. For this case, we set the 
maximum CFL number to 50, and were able to gain the cost reduction factor of 0.127 com- 
paring with fully implicit algorithm. Recall that the cost reduction factor was defined in 
Section 5 as the ratio of the cost marching in time (and converging to steady state) between 
the implicit/explicit and fully implicit algorithm. Thus, the above factor indicates that the 
implicit /explicit algorithm was about 8 times faster— a considerable speedup. 
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Figure 7.57: Blunt body problem, initial conditions and geometry. 
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Figure 7.58: Blunt body problem, initial computational mesh. 
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Figure 7.60: Blunt body problem, density contours for the inviscid solution. 
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Figure 7.61: Bl un t body problem, final mesh with additional manual adaptation in the near 
wall region. 
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Figure 7.62: Blunt body problem, density contours for the final mesh. 
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Figure 7.64: Graded Mesh, includes 9 layers of quadratic elements along the cylinder (Re- 
ferred as Mesh B.) 
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Figure 7.65: Mesh B after 2 levels of adaptation. Elements on the wall are of the second 
order. 
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Figure 7.66: Pressure contours obtained on the adapted Mesh B. 
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Figure 7.67: Sonic line obt&ined on &diptcd Mesh B. 
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Figure 7.68: 3D view of density in the shock impingement region. 
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Figure 7.69: 3D view of pressure in the shock impingement region. 
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Example 10: Shock /Boundary Layer Interaction With Separated Flow 

The second benchmark problem involves a shock /boundary layer interaction with separated 
flow (Holden problem [30,32]). The particular flow conditions used in the analysis of the 
Holden problem are as follows: 

M = 14.1 

Re = 72,000 

= 80 °K 

T wan = 297.2 °K 

inclination of the wedge = 24° 

We initiated work on this problem in the previous phase of the project, but encountered 
considerable difficulties in resolving the flow recirculation region. These difficulties were 
resolved in the Phase II effort and, moreover, the source of our previous difficulties was 
identified as the artificial dissipation model. The effort spent on solving this problem during 
the previous phase is summarized in the next two paragraphs, and the related numerical 
results are included in Figs. 7.72 to 7.77. 

The first approach for solving this problem was with the two-step algorithm using an 
initial mesh of 11 x 25 linear elements. The elements were stretched in the horizontal 
direction with the aspect ratio along the solid wall boundary of only about 1:5 and the size 
of the smallest elements is about 0.03. The problem turned out to be practically unsolvable 
without the artificial viscosity especially designed for highly stretched or distorted elements. 
The standard artificial viscosities (Lapidus’ and Morgan’s models described in Section 2.4) 
failed to stabilize the solution around the stagnation point. On the other hand, introducing 
elements with an aspect ratio of 1 in the wall area results in a prohibitively large number 
of elements and small time step, slowing down the integration process. However with our 
modified artificial dissipation model the solution process was able to proceed on the original 
mesh. An /i-adaptive mesh with three levels of refinement (based on residual error indicators) 
is shown in Fig. 7.72, the corresponding density contours are shown in Fig. 7.73. Salient 
points about the solution to note include: sharp shock resolution, minimal number of degrees 
of freedom to capture the shocks, and reflected shock/leading edge— shock/boundary layer 
interaction. However, the key viscous feature along the wall, the recirculation bubble, has 
not been resolved. 

In a parallel modeling effort, we also used the one-step Taylor-Galerkin algorithm which 
allows one to use much larger time steps and therefore leads to faster convergence of the 
steady state solution. The initial linear mesh used in this solution procedure is shown in 
Fig. 7.74. Compared to the previous mesh we introduced a considerably finer discretization 
along the solid wall boundary, the size of the smallest elements is now 0.005 and the maximum 
aspect ratio is 25. The problem was run with the CFL constant set to 2.5. From the solution 
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obtained after 100 and 200 time steps, the mesh was /i-adapted up to the second level, based 
on the residual error indicators, see Fig. 7.75. The contours of the density for this mesh and 
v-component of the velocity are shown in Figs. 7.76 and 7.77. While the mesh also provides 
a reasonable resolution of the shocks, the recirculation region was still not well developed. 
For this mesh we also introduced a p-enriched mesh along the solid wall, with the hope for 
better resolution of recirculation of the flow, but this effort was also unsuccessful. 

During phase II of this project, we have experimented with several different artificial 
dissipation (AD) models. We found that this is the crucial factor in resolving the flow 
separation phenomenon for this benchmark problem. Recall that in Section 2.4 we have 
implemented three different AD models: 1) the classical Lapidus, 2) the modified Lapidus 
by Lohner at el., and 3) our version which is a modification of the second model in order 
to handle elements of high aspect ratios. As expected, out modified model performs better 
than the other two models on elements with irregular shapes (high aspect ratios or high 
distortions), and we have used this model successfully to solve all the test problems, in- 
cluding a compression corner problem with a lower Mach number as described previously. 
For Holden’s problem at Mach 14, however, all three models have difficulty resolving the 
separation region. Usually, in the solution of problems with strong shock-boundary layer 
interactions, not only does the mesh have to be well adapted according to the flow features, 
but the AD is also a key factor. From numerous test runs, we have found that for this 
problem, resolution of the recirculation bubble is extremely sensitive to the application of 
AD because of the flat shock angle and sharp boundary layer pattern involved. It is well 
known that, for viscous compressible hypersonic flows, the ideal AD model should provide 
just enough dissipation to smooth shock discontinuities without oversmearing, and should be 
formulated to avoid introducing too much dissipation (relative to the natural dissipation) in 
the boundary layer region. This means that the AD should have an accurate built-in sensor 
to control the amount of dissipation and detect the direction in which the AD is to be added. 
All three models (belonging to the Lapidus family) are using the velocity gradient as sensors, 
and their success are based on the argument that the maximum change of velocity gradient 
is perpendicular to shocks so the AD is also added in that direction. In the boundary layer 
region, since the change of velocity gradient is perpendicular to the solid wall, there will be 
no AD introduced in the tangential direction which may otherwise put too much “artificial 
dissipation to accurately resolve the “real” flow features. 

In the formulation of our AD model described in Section 2.4, the unit vector t was 
based on the gradient calculated on the master element as in Eq. (2.95). Noticing that the 
orthogonality is not preserved under the transformation, we used the l vector computed 
on the original element. With this simple remedy, the recirculation bubble was successfully 
resolved. However, because it is a less dissipative mechanism than other models, the solution 
is more oscillatory near the shocks and the leading edge, and smaller time steps were required 
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to successfully converge. The analysis effort based on these various AD models is presented 
below. 

In order to compare the numerical results with those presented in ref [30], we used an 
initial mesh consists of 27 by 25 linear elements, see Fig. 7.78, which was clustered almost 
identically to the SM mesh used therein. After two levels of uniform refinement, the thickness 
of the smallest elements would be 2.42 * 10 ~*ft, which also matches the SM mesh. On the 
mesh with one level of uniform refinement, we have experimented with three different AD 
models: (a) the Lapidus, (b) the modified version of Morgan’s presented in Sec. 2.4, and 
(c) our modified model with the fix as described in the previous paragraph. The results are 
shown in Figs. 7.79, 7.80, and 7.81, respectively, and are displayed by the contours of the v- 
component of the velocity. While the first two cases failed to resolve the recirculation bubble, 
the last one clearly indicates the flow separation region. At this point it is important to note, 
that we do not claim here that this variant of artificial dissipation is ultimately better than 
others. It rather seems that the other models were too dissipative for this problem and were 
“wiping out” the fragile separation point. The design of an ultimate artificial dissipation 
model, which would resolve the shocks without smearing other features of the solution, is 
still an open challenge in computational fluid dynamics. 

The same mesh with two levels of ^-adaptation is shown in Fig. 7.82. The numerical 
results presented in Figs. 7.83 to 7.85 show contour plots of the density, v-component of 
the velocity, and a 3D view of pressure in the recirculation region. A comparison of the 
pressure, skin friction, and heat transfer coefficients along the solid wall with the experiment 
data are profiled in Figs. 7.86 to 7.88, respectively. In general, they show a similar overall 
agreement with experiment data as those numerical results for mesh SM presented in ref [30], 
except that our heat flux is underpredicted in the shock reattachment region. We believe 
that this discrepancy is caused by the AD introduced in this region. Intuitively speaking, a 
certain fraction of the total dissipation on the wall is “taken over” by the AD model, which 
tends to reduce the flux contributions from the natural dissipation. Note that, since the 
amount of dissipation due to the AD model decreases with decreasing mesh size, further 
mesh refinement would produce even better agreement between the numerical results and 
experimental data. 

In the last stage of these computations, the recently implemented implicit/explicit al- 
gorithm was used for the solution of the problem. Due to the oscillatory behavior of the 
solution (caused by the high speed flow) near the leading edge, the maximum bound for 
CFL number was set to a rather low value of 2. On the final adapted mesh this corre- 
sponds to only 6.6% of domain selected as implicit, with the implicit elements clustered in 
the near-wall region. Based on this zone selection, the cost reduction factor with respect 
to the fully implicit and fully explicit schemes was 0.06 and 0.65, respectively. This means 
that the implicit /explicit algorithm converged about 17 times faster than the fully implicit 
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algorithm (the only method available previously in the code) and about 2 times faster than 
the fully explicit algorithm, with am additional beneficial effect of smoothing the oscillatory 
tendencies of the solution near the plate tip. 
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Figure 7.72: Holden problem, Re = 72,000, h-adapted mesh with three levels of refinement. 
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Figure 7.73: Holden problem, Re = 72,000, density contours for an /i-adapted mesh. 
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Figure 7.74: Holden problem, Re = 72,000, initial mesh used for the one-step method. 
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Figure 7.75: Holden problem, Re = 72,000, A-adapted mesh. 
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Figure 7.76: Holden problem, Re = 72,000, density contours for an A-adapted mesh. 
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Figure 7.78: Initial mesh (referred as SM mesh). 
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Figure 7.79: V-velocity contours for Lapidus AD model (on SM mesh with one level of 
uniform /i-refinements. 
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Figure 7.81: V*velocity contours for modified 
with one level of uniform A-refinement. 
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Figure 7.82: SM mesh after 2 levels of /(-adaptation. 
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Figure 7.85: 3D view of pressure distribution for an adapted SM mesh (the protion of the 
mesh displayed covers the separation and reattachment regions. 
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Example 11: Rearward Facing Step with Strong Expansion 

The third benchmark problem modeled herein is a supersonic, viscous, laminar separated 
flow over a rearward facing step. A schematic of the geometry is shown in Fig. 7.89 with 
several additional labels identifying regions of the domain which are referred to later. The 
experimental data provided by NASA Langley for this test case, as well as several other 
similar configurations, is found in Ref. [17]. In this reference, experiments were conducted for 
a range of test conditions: T , ^ = 2700 ~ 5500°#, Me© = 3.95 ~ 4.27, Rego/cm — 160 ~ 2200. 
For the non-suction case, with no chemical reactions which is our focal point in this example, 
the experimental data includes a heat transfer distribution and a surface pressure distribution 
downstream of the step. 

The experimental setup contained two sets of control data from which we selected the 
case with the maximum step height. For this case, h = 1.02 cm, and the following flow 

conditions exist: 


Moo = 4.08 

Rtoo,h - 1650 

Too = 3750° K 

r waU = 297° A' 

The initial mesh for this problem consists of 28 by 8 linear elements with some clustering 
near the wall as shown in Fig. 7.89. The initial conditions were defined as a uniform flow 
in the whole domain, and the problem was first solved with 2 levels of uniform refinements 
which corresponds to 3729 degrees of freedom. A very small time step was initially required 
to avoid the negative densities and pressures near the corner of the backstep during start up 
from free stream conditions. After the flow stabilized, we applied the ^-adaptation option to 
begin resolving the shocks and boundary layer regions. The flow features of specific interest 
in this problem required high resolution of the flow field variables in the regions of the 
boundary layer along the wall, and around the two corners of the backstep. 

The corresponding mesh with 3 levels of A-refinement is shown in Fig. 7.90, and a 3D 
view of the density is presented in Fig. 7.91. Large gradients of the solution occur near 
the wall both upstream and downstream of the backstep. For the upstream portion of the 
computational domain there are approximately 5 layers of Unear elements in the boundary 
layer which are probably adequate to roughly resolve the flow features therein. In the 
downstream portion of the computational mesh, however, there are only approximately two 
layers of Unear elements which are undoubtedly insufficient to capture the heating rates and 
other fine scale phenomena. 

To adequately resolve the downstream portion of the mesh we applied a uniform h- 
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refinement to four layres of elements near the wall in this region to arrive at the mesh 
(referred to as M3-mesh) as shown in Fig. 7.92. On this mesh we continued the solution 
process with a time-step size corresponding approximately to CFL = 5. The resulting 
pressure contours on the entire mesh and the temperature contours in the backstep and 
downstream regions are shown shown in Figs. 7.93 and 7.94, respectively. From these plots 
one can observe the overall patterns of the leading edge shock, the expansion fan at the corner 
of the step, the reattachment shock downstream of the step, and the boundary layers near 
the wall. The density and heat flux coefficient along the wall downstream of the backstep 
are also profiled in Figs. 7.95 and 7.96, respectively. 

From this point we continued the solution process along two different but parallel paths 
to investigate the effectiveness of fi-adaptation versus p-enrichment procedures. The first 
approach consisted of applying the /i-refinement /unrefinement algorithm to the M3-mesh 
using the element residual error estimate scheme presented in Sec. 4.1 as the adaptation 
driver. This adaptive pass resulted in the refinement of several elements in the following 
regions: the leading edge, the boundary layers near the wall upstream and downstream of 
the step, the corner of the step, and the reattachment shock zone. The resulting /i-adapted 
mesh (referred to as M4H-mesh) which consists of 9700 linear elements and 9633 degrees of 
freedom is shown in Fig. 7.97. After continuing to march the solution for approximately 6 
seconds of real time, the solution reached steady-state on the new mesh. The corresponding 
pressure contours, temperature contours, density profile, and heat flux coefficient are shown 
in Figs. 7.98 to 7.101. As expected, the enhanced mesh provides a better resolution of the 
flow variables in the reattachment shock and boundary layer regions than those obtained on 
the M3-mesh. 

The second approach consisted of using the p-enrichment procedure proposed in Sec. 4.2 
to adapt the mesh. Based on the solution obtained on the M3-mesh, the mesh was first 
adaptively /i- unrefined (as in the M4H-mesh) to remove excess degrees of freedom. Then, 
for the elements with errors larger than the user specified threshold value (the same as the one 
used for the M4H-mesh), we applied the adaptive directional p-enrichment procedure. The 
directional adaptation indicators were calculated based on the derivatives of the solution, 
Eq. (4 42), and the values of and k for selecting the enrichment directions were set to 0.3 
and 0.7, respectively. The enriched mesh (referred to as M4P-mesh) is shown in Fig. 7.102. 
Figs. 7.103 and 7.104 show blowups of the mesh around the corner of the step and along 
the wall downstream of the step, respectively. Notice that the M4P-mesh consists of only 
7774 degrees of freedom, which is almost 2000 less than the M4H-mesh. After marching 
the solution for 6 seconds of real time on this mesh, the corresponding pressure contours, 
temperature contours, density profile, and heat flux coefficient were extracted as shown in 

Figs. 7.105 to 7.108. 

Qualitatively, the p-enriched M4P-mesh showed a similar overall improvement in the 
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solution over M3-mesh as the /i-refined M4H-mesh did. The density profiles along the wall 
downstream the step of M4H-mesh and M4P-mesh appear to be virtually the same. The one 
significant difference in these two results appears in the prediction of heat flux coefficients 
which shows about a 20 to 30 percent variation. Since both cases predicted the same location 
of maximum heat transfer, we believe that the descripency of wall heat transfer rate in the h- 
adapted mesh is caused by the still insufficient mesh clustering in the direction normal to the 
boundary layer near the wall. Based on our experience in solving the other two benchmark 
problems (the Holden’s problem and the blunt body with impinging shock problem), where 
the heat transfer rates were always underpredicted on fc-meshes, we tend to conclude that 
the p-enriched mesh in the boundary layer region is providing a more effective use of degrees 
of freedom than /i-refined mesh. 

It should be noted that for this benchmark problem the experimental data presented in 
Ref. [17] do not provide enough detailed information to compare our numerical results with. 
Furthermore, unlike the other reference (e.g. [19]) where real gas effects were taken into 
account, our assumption of perfect gas at 7 = 1.4 may also make the comparison between 
numerical and experiment data difficult at this time. Note thar the real gas effect will be 
considered in the next year of this project. 
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Figure 7.90: Adapted mesh with 3 levels of h-refinement. 
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Figure 7.91: 3D view of < 
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Figure 7.92: Adapted mesh with uniform /i-refinement near the wall downstream of the 
backstep (M3- mesh). 
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Figure 7.96: Heat flux coefficient profile along the wall downstream of the backstep, 
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Figure 7.97: The final A-adapted mesh (M4H-mesh). 




Figure 7.9S: Pressure contours, M4H-mesh. 
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Figure 7.99: Temperature contours in the backstep region, M4H-mesh. 
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Figure 7.100: Density profile along the wall downstream of the backstep, M4H-mesh. 
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Figure 7.101: Heat flux coefficient profile along the wall downstream of the backstep, 
M4H-mesh. 
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Figure 7.102: The final p-adapted mesh (M4P-mesh). 
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Figure 7.103: Blowup of M4P-mesh around the top corner of the backstep. 
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Figure 7.104: Blowup of M4P-mesh near the wall downstream of backstep. 
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Figure 7.106: Temperature contours in the backstep region, M4P-mesh. 
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Example 12: Double Swept Wedge Comer Flow Problem, M = 3 

During the last year of the project, we have also initiated the solution of a three-dimensional 
benchmark problem, which involves modeling of the invicid flow past a wedge consisting 
of two planes inclined at different angles to the direction of the flow. The geometry of 
the problem is rather simple and is shown in Fig. 7.109. Due to plannar symmetry of 
this geometry the problem was solved only in one half of the computational domain with 
appropriate symmetry boundary conditions. 

Up to date, we have only performed an initial solution of this problem on a rather coarse 
mesh. The corresponding solution, contours of density, and an ^-adaptive mesh of linear 
elements are presented in Fig. 7.110. The flow, as expected is characterized by a skew shock 
which leaves the computational domain without any reflexions. Obviously, the results are 
far from converged on such a coarse mesh, however, the general shock structure and other 
features of the flow appear to be developing correctly. Further computations for this and 
other three-dimensional problems will be performed in the next year of the project, after 
completion of development of efficient three-dimensional implicit/explicit algorithms. 
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Figure 7.110: Double swept wedge, M = 3, density contours and A-adapted mesh of linear 
elements. 
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8 Phase II Project Summary and Future Directions 

The computational results obtained for the various test cases and the theoretical advances 
made in the area of adaptive finite element methods over the past four years has been 
quite encouraging. They indicate that the h-p finite element method is not only a feasible 
approach for solving hypersonic flows but when combined with an implicit/explicit solution 
methodology provides an optimal framework for systematically changing the structure of the 
computational mesh to provide highly accurate numerical results with a minimal number of 
degrees of freedom and at minimum computational cost. 

The technical efforts over the past year have focused on two topics which are specifically 
related to the overall performance of the flow solver. The first of these issues is that of 
directionally-dependent error estimation schemes which in general focused on an automated 
procedure for choosing an appropriate direction for p-enrichment. The algorithm selected 
herein uses a residual error estimate in conjunction with gradients of the residual error esti- 
mator or gradients of the solution (either of which may be selected by the user.) The residual 
error estimate itself is used to identify elements of the computational domain with relatively 
high errors. Among the group of elements with high errors a directional pointer, based on 
the gradients of a specified quantity, is obtained which indicated an optimal direction (or di- 
rections) in the master element for p-enrichment. Several numerical results have shown that 
the algorithm in general selects directions for enrichment that are normal to boundary layers 
and/or normal to shocks. In regions where point types of singularities exist, the indicator 
tends to select isotropic types of refinement. 

The second major topic addressed over the past year, which is related to the overall 
performance of the flow solver, is that of implicit/explicit solution procedures. The general 
idea behind this methodology is that there is often a large diversity of stable time steps for 
various parts of the computational mesh. If one uses an implicit method in regions of the 
mesh with relatively severe time step restrictions and an explicit solution method in other 
regions of the mesh where the stability restrictions are not as severs then an optimal balance 
of computational effort is achieved within a single time step. Such an algorithm based on 
a general class of implicit Taylor- Galerkin algorithms was implemented within the context 
of the two-dimensional h-p flow solver. Based on the results of several test problems which 
employed this algorithm, an average computational savings of about 25 percent was achieved 
when compared with fully explicit algorithm and about 60 percent when compared with fully 
implicit algorithm used previously. Note that higher computational savings were obtained 
on problems with larger variations in the mesh size (see Section 7 for specific details). 

The results obtained over the past year for the test problems and the benchmark problems 
in general indicate that the implicit/explicit methodology is a key component of an efficient 
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solution process. It provides a second level of optimization whereby not only an optima] mesh 
is used in the solution sequence but also an optimal time stepping procedure is employed. 
The proposed next phase of the effort will focus on extending the implicit/explicit solver 
methodology and the directional error estimation strategies to the three dimensional case. 
Based on our results over the past year the effort should be highly successful and may even 
provide more computational savings than in the two dimensional case. 
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Appendix 

A Performance Issues 

This appendix contains a summary of the work performed at the Computational Mechanics 
Company, Inc. in the areas of iterative solution techniques for the h-p finite element method 
and static condensation for high order spectral elements. The goal of these studies was to 
improve the computational efficiency of the h-p finite element method through improved 
element conditioning and reduction of element storage requirements. Funding for this work 
was provided through NASA, the Department of Defense, and private funding, and was not 
a part of the research and development effort conducted as part of this project. The results, 
however, are related to this project and are included here for completeness. 

A.l Iterative Methods for the h-p Finite Element Method 

This section presents some numerical studies on the performance of iterative solvers imple- 
mented in the h-p finite element environment. In particular, a study of the GMRES algorithm 
combined with various versions of the block Jacobi preconditioners was performed. This iter- 
ative solver is used to solve nonsymmetric systems of linear equations occurring in the finite 
element analysis of compressible and incompressible flows, as well as in certain problems in 
solid mechanics. 

The most difficult issues in the iterative solution of these systems is related to the defini- 
tion and implementation of the preconditioner. A good preconditioner must satisfy several 
somewhat contradictory requirements, in particular: 

• It should improve the conditioning of the system to be solved. 

• It must precondition out all the penalty terms used to enforce the boundary conditions 
and other constraints. This is because extremely large eigenvalues introduced by the 
penalty terms ruin the convergence of the GMRES accelerator. 

• The submatrices inverted in the process of preconditioning should be relatively small 
to minimize the computational cost of preconditioning. 

In order to reasonably compromise the above criteria, we have studied a few different pre- 
conditioners of the block Jacobi type. These preconditioners were based on the concept of a 
“patch,” which is explained further in this section. 

The following subsections discuss the basic algorithm, preconditioning, and numerical 
studies on the performance of the GMRES method. 
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A.1.1 Basic Description of the Iterative Method 

GMRES is a method for solving non-symmetric linear systems of equations, as described in 
[23]. The following is a brief description of the method: 

Consider the system of linear equations: 

Ax — f 

The algorithm of the truncated/restarted GMRES method consists of the following steps: 


1. Start: 

Choose the initial guess xo 
Compute the search direction: 

r 0 = / — Ax o 


normalize: 

v,= ra 

2. Loop through the Gramm-Schmidt algorithm to calculate m Arnoldi vectors spanning 
the Krylov subspace. 

This is done by the iterative procedure: 


j 

= 1,2 ... m 

h-ij 

= ( Avj,vi ) ; * = 1,2 — j 

V)+l 

j 

= A.Vj ^ 
«=1 

^j+1 0 

= IIVj+. II 

Vj+l 

^J+lJ 


3. Form the approximate solution 

X m — Xq 4’ ^m!/m 

where y m minimizes ||& - y € iT (see [23] for details). 
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4. Restart: 


Compute r m = / — Ax m ; if the residual is sufficiently small, then stop 
else compute: 


x 0 = 

Vi = 

and go to 2. 

A. 1.2 Separating a Preconditioner From an Accelerator 

The GMRES algorithm is used to solve a preconditioned system of equations: 

D-\Kx-f) = 0 

The idea of using a preconditioner is to replace the original system of equations Kx — f = 0 
by the new one, as above, with the matrix D~ l K characterized by a much lower condition 
number than that of K. 

The preconditioner D can be constructed from any so-called “basic iterative method” as 
follows: 

Given a basic iterative method for solving Kx — f — 0: 



In+l Gin + k 


where G, k satisfy the consistency condition: (G — /)x + k = 0 is equivalent to Kx = /. We 
introduce D such that the above formula becomes: 

x n+1 =D-'[-(K-D)x n + dk] 

Therefore G = —D~ 1 (K — £>)—► D~ l = (I — G)K~ 1 . In many situations, even when the 
basic iterative method does not converge, it can serve as a good preconditioner (such as the 
. Jacobi algorithm applied to the mass matrix in three dimensions). 

Solving the preconditioned system of equations D~ l (Kx — /) = 0 with an algorithm like 
GMRES requires performing the following operations: 

y : = D-'Kx 

r : = D-'(-Kx + f) 
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From the definition of D we find that instead of explicitly preconditioning the matrix K by 
D'\ one can perform one step of the basic method to compute: 

y = x n - x„+i with / = 0 
r = x„ — x„+i with the actual / 


This means that the “basic iterative method” = Gx„ + k can be quite independent from 

the accelerator GMRES and the above operations can be performed just by c mg rom 
GMRES a “basic iterative solver,” block Jacobi or Gauss-Seidel in our case. 


A.l.S Patch Definition 

The solver in particular the block Jacobi preconditioner, uses a block of stiffness matrices 
which are associated with patches of elements. A patch is the assembly of nodes associated 
with a node being preconditioned. The block inverted in D~ includes all the degrees 
freedom nodes included in the patch. Three definitions of patches have been investigated: 


• patches associated with global linear nodes, 

• patches associated by single nodes, and 

• a combination of (i) corner nodes on the boundary and (ii) interior nodes. 


The following is a brief description of each of the patch definitions: 

1 Patches associated with global linear nodes consist of all nodes connected with the 
central node through the edges of elements. An example of a patch 1 is presented in 

Fig. A.l, where: 

• = patch node 
x = activated node 

2. Patches associated with single nodes (Fig. A.2). An active node is the patch node. 

3 Version 3 combines patch 1 for boundary nodes and patch 2 for interior nodes. An 
' example of this type of patch is presented in Fig. A.3 where . = patch node and x - 

activated node. 


Patches of the first kind provide a quite powerful preconditioning. Its major drawback 
is the large size of the patch stiffness matrix, especially for high p and thre^dimensional 


elements. 
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• = patch node 
j = activated node 
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Figure A.l: Patch type 1. 



Figure A. 2: Patch type 2. 




Figure A. 3: Patch type 3. (a) Boundary Node, (b) Interior Node. 
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Patches of the second kind introduce small stiffness matrices. It is, however, not very ef- 
fective, as it does not precondition out the penalty terms used to enforce boundary conditions 
for higher orders of approximation. 

Patches of the third kind combine smaller patches on the interior and larger patches on 
the corner nodes of the boundary. It does precondition against the penalty terms on the 
boundary, however, it is not as effective as patches of the first kind. Note that for linear 
elements, patches of the first, second, and third kind all produce the same results. 

A. 1.4 Numerical Results 

In this section, some sample results are presented from our study on the performance of 
iterative methods (i.e., GMRES accelerator and block Jacobi preconditioner). This perfor- 
mance study consists of a convergence check and a timing study where a comparison with a 
direct frontal solver is made and also the effects of using different kinds of patches is studied. 
The test problems considered range in complexity from a simple Lj projection to problems 
in fluid mechanics (both incompressible and compressible). 

Example 1. L 2 Projection 

This problem is a two-equation I 2 projection case: x 2 + y 2 = 1.0 on a square domain. It 
is the simplest of the test cases and is used to verify the basic performance of the iterative 
solver. The mesh consists of second order polynomial element shape functions. There are 4 
elements and 25 degrees of freedom. An example of the mesh is shown in Fig. A.4. 

From Table A.l it can be seen that only the iterative solver using patches of the first 
kind converges. It takes 2.5 times longer to converge as compared to the frontal solver. This 
is not a surprising result, however, because for small patches, direct solvers usually perform 
better than iterative solvers. 

Example 2. Incompressible Fluid Mechanics 

Here we present the timing results for a two-dimensional driven cavity flow. The geometry 
and boundary conditions of the problem are simple and are very well known. For the case 
under consideration, the flow field moves in the positive x-direction and the Reynolds number 
is 1.0. The mesh consists of second order polynomial element shape functions for the velocity 
field and first order polynomial element shape functions for the pressure field. There are 64 
elements and 289 degrees of freedom per unknown. An example of the mesh is presented in 

Fig. A.5. 
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Figure A. 4: Mesh for the Xj projection. 
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Table A.l: Timing results for Xj-projection. 
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Figure A. 5: Driven cavity mesh. 
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Table 2 

Timing results for driven cavity problem. 
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From Table A. 2 it can be seen again that only the iterative solver with patches of the 
first kind converges. It takes approximately 2.9 times longer to reach the same solution as 
the frontal solver. The iterative solver with patches of the second and third kind do not 
converge to the correct solution; however, the residuals of the GMRES does converge. As 
expected, the iterative solver with patches of the third kind converge closer to the solution 
as compared to patches of the second kind. The bad performance of preconditioners based 
on patch types 2 and 3 is probably caused by a complex coupling of pressures and velocities 
(through the gradient operator), which is not resolved by small patches for interior nodes. 

Example 3. Compressible Fluid Mechanics 

Presented here are the timing results for the standard wedge benchmark problem for inviscid 
supersonic flows. For this case, a 5 degree wedge with a uniform inflow and the Mach number 
of 5.0 was selected. The mesh consisted of 240 linear elements with 275 degrees of freedom 
per unknown. (The mesh is presented in Fig. A. 6.) The problem was run for 300 time steps. 

From Table A. 3 it can be seen that there is no convergence check. This is due to the 
frontal solver and the iterative solver being converged to the same solution. As expected, 
the iterative solver converged irrespective of the type of patch used. The iterative solver 
using patches of the first kind converged the fastest followed by patches of the second kind, 
and lastly, patches of the third kind. The iterative solver converged 3.9 times faster than 
the frontal solver. The primary reason the iterative solver converged faster than the fronter 
solver is that the mass matrix for this case is symmetric and positive definite. A symmetric 
positive definite matrix has real eigenvalues and has a good conditioning number. These 
characteristics of the matrix bring about a faster convergence rates for iterative solvers. 

Conclusions 

The iterative solver (i.e., GMRES accelerator with block Jacobi preconditioner), utilizing 
patches of the first kind, converged to the same solution as the frontal solver for all cases 
under consideration. The iterative solver converges slowly if the matrix is ill-conditioned or 
if a penalty method applied to the boundary condition has not been preconditioned out. It 
was found that the iterative solver converges more rapidly than the frontal solver for well- 
conditioned matrices (the presence of mass matrix contributions on the left hand side is very 
beneficial). From the present study, the use of the iterative solver with patches of the second 
and third kind does not appear to be a good option. 
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Figure A. 6: 5 degree wedge mesh. 
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Table A. 3: Timing results for 5 degree wedge problem. 
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A.2 Matrix Condensation 


One of the side effects of using higher order finite element approximations is the extremely 
large size of the element matrices that result, especially in three dimensions. For an eighth 
order element in three dimensions, this corresponds to a local vector with 729 x NDOF where 
NDOF is the number of solution components. Converting this into storage requirements 
implies that approximately 0.5 megawords of memory axe required for each component of the 
solution for a single element. To reduce this large memory requirement we are investigating 
the use of static condensation to eliminate the internal degrees of freedom before the element 
matrix is passed to the iterative solver. For the eighth order element mentioned above, this 
corresponds to the removal of 343 x NDOF degrees of freedom at the element level which 
in turn reduces storage requirements of the element matrix to less than one third of the 
original size. In addition to reducing the element storage size, this technique also has the 
advantage of improving the conditioning of the global system by eliminating the coupling 
between nodal and central degrees of freedom and edge and central degrees of freedom. 


Condensation Algorithm 


Depending on which numerical algorithm is actually used in the solution of the governing 
equations, the resulting forms of the element local stiffness matrix and right hand sides may 
be written as 


Si = 


Ai B, 

Di Ci 


R, = 


hi 


Ci 


where 

A 


Ci 


B x and D, 


bi 

Ci 


are stiffness matrix components representing interactions among the 
internal degrees of freedom only. 

are stiffness matrix components representing interactions among the 
side and nodal degrees of freedom. 

are stiffness matrix components representing mixed types of interac- 
tions between internal and side and nodal degrees of freedom. 

is the right hand side vector associated with the internal degrees of 
freedom. 

is the right hand side vector correspnding to all other degrees of freedom 
not at the central node. 


The process of condensation, 
Schurr complement: 


or elimination of internal unknowns, entails computing the 

Ci = Ci - D x A~'Bi 
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and modifying the right hand side as follows 

Ci = Ci- DiA-'bi 

These quantities are then normalized so that all local diagonal matrix entries are one. 

Ci = E~^C{E~i Ci-E~^Ci where Ei = diag (c.) 

Using these element matrices, a global system of the following form is obtained which no 
longer contains the internal degrees of freedom, 

Sa = r 

where 

S = Y J C i r = '$2 Ci 

i i 

and 5i is a vector of unknowns associated with the side and nodal degrees of freedom. 

The final step is the calculation of the internal unknowns a;. These are computed using 
5, by solving the local systems 

AiCti — bi — Bidti 

Matrix Operations and Array Sizes 

To solve a problem of q equations on a n x n mesh using uniform p order elements, the size 
of the global stiffness matrix is: 

[((n + l) 2 + 2n(n + l)(p - 1) + n 2 (p - l) 2 ) x 9] 

where (n + l) 2 represents the degrees of freedom associated with the combinations of the 
corners, 2n(n + l)(p — 1) represents the degrees of freedom associated with the combinations 
of the edges, n 2 (p — l) 2 represents the degrees of freedom associated with the center, and q 
is the number of equations. 

For an 8 x 8 mesh with fourth order elements and two equations, one solves a system of 
2178 equations. With condensation, this system is reduced to n 2 inversions of [(p — l) 2 x q ] 2 
matrices and to solving a system of [(n + l) 2 + 2n(n + l)(p — 1)] x q equations. For the 
sample case, this is 64 (18 x 18) matrices to invert which is equivalent to solving a system 
of 513 equations. Locally it reduces the stiffness matrices from (p+ l) 2 x q to 4 pq, which 
is significant for high p. It is even more significant in three dimensions when it goes from 
(p + l) 3 x q to (6p 2 + 2 )q. 

Preliminary results using the static condensation procedure for L 2 projections has shown 
a slight increase in the converge rate which is most probably due to improved conditioning 
and up to 20 percent savings in computational time for larger p. 
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